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Abstract 


In  a  recently  developed  material  forming  method  called  Functionally  Gradient 
Materials  as  well  as  applications  of  material  deposition  processes  such  as  ion  plating, 
composite  materials  are  being  created  where  the  interface  possesses  a  gradually 
varying  material  composition  and  properties.  This  study  was  directed  at  the 
mechanics  of  such  materials  when  there  is  a  crack  on  the  interface,  and  sought  to  find 
parameters  that  govern  the  crack  growth  such  as  the  crack  tip  stress  intensity  factors, 
strain  energy  release  rate  and  probable  direction  of  crack  extension.  The  mixed 
boundary  value  problem  involved  two  bonded  materials  having  finite  thicknesses  with 
an  interface  crack  under  plane  strain  or  generalized  plane  stress  conditions.  One 
material  is  homogeneous  and  the  other  nonhomogeneous  with  an  exponential  property 
variation  in  the  y -direction. 

Fourier  transforms  was  applied  to  Navier’s  equations  to  derive  a  system  of 
singular  integral  equations  with  a  simple  Cauchy  kernel  and  Fredholm  kernels.  The 
x -derivatives  of  the  two  crack  opening  displacements  are  assumed  to  be  the  unknowns. 
Extensive  asymptotic  expansions  of  the  kernels,  which  were  the  algebraic  sum  of 
rational  functions  of  7  by  7  and  8  by  8  determinants  as  the  numerator  and 
denominator  were  carried  out  in  order  to  separate  the  Cauchy  kernel  and  to  facilitate 
the  integral  equations  numerical  computation.  The  problem  was  solved  numerically 
by  converting  to  a  system  of  linear  algebraic  equations  and  by  using  a  collocation 
technique.  The  stress  field  near  the  crack  tip  is  mixed  mode  and  is  shown  to  have  a 
standard  square  root  singularity. 

Cases  for  a  wide  range  of  degrees  of  nonhomogeneity  and  combinations  of 
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material  thickness  to  crack  length  ratios  were  computed  with  loadings  of  uniform 
normal  stress  and  uniform  shear.  The  results  of  a  special  case  where  both  materials 
have  very  large  thicknesses  compare  very  well  with  the  previous  results  for  two 
bonded  half  planes. 

The  technique  developed  in  this  study  is  useful  for  fracture  mechanics  studies 
of  composite  materials  that  have  a  nonhomogeneous  interfacial  zone.  The  results 
compiled  are  especially  well  suited  for  studying  the  delamination  problem  in 
nonhomogeneous  thin  films  bonded  to  elastic  substrates. 
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Chapter  1 
Introduction 


1.1  Introduction 

Materials  of  composite  nature  have  found  their  applications  perhaps  as  early 
as  the  recorded  history.  One  of  the  earliest  applications  may  have  been  in  protection 
against  corrosion  by  using  a  variety  of  coatings.  The  majority  of  the  more  recent 
applications  of  composite  materials  seek  to  utilize  the  strengths  or  mechanical 
properties  of  component  materials  in  order  to  optimize  the  properties  of  the  combined 
system.  For  example,  blades  in  gas  turbine  engines  are  subject  to  high  stresses  and 
highly  corrosive  environment  of  oxygen,  sulfur  and  chlorine  gases.  A  monolithic 
material  such  as  a  high  temperature  alloy  is  incapable  of  providing  both  functions. 
The  solution  is  to  design  the  bulk  for  the  mechanical  properties  and  use  a  coating  of 
another  alloy  to  provide  the  corrosion  resistance.  Another  example  of  composite 
materials  is  dry  film  lubricant  coatings,  as  in  the  case  of,  for  example,  using  thin  gold 
films  bearing  and  sliding  parts.  Such  dry  film  lubricants  are  important  for  critical 
parts  where  conventional  organic  fluid  lubricants  are  susceptible  to  degradation. 

In  the  quest  of  finding  an  optimal  composite  for  a  particular  need,  the  problem 
of  thermal  mismatch  at  the  bimaterial  interface  has  always  been  one  of  the  major 
reasons  that  cause  the  material  to  fail  to  perform  its  intended  function.  Thermal 
mismatch  and  residual  stresses  occur  because  the  temperature  under  which  the 
composites  are  processed  does  not  fall  in  the  temperature  range  where  it  is  supposed 
to  perform,  and  usually  far  from  it. 

With  the  help  of  novel  material  forming  processes  and  innovative  apparatuses, 
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material  scientists  have  found  ways  to  produce  new  materials  with  maximum 
efficiency  as  never  before.  In  laboratories  and  manufacturing  facilities,  physical  vapor 
deposition  (PVD),  chemical  vapor  deposition  (CVD)[1]  and  their  derivatives  have 
been  used  to  create  materials  from  the  very  elemental  building  blocks,  i.e.  atoms  and 
molecules  for  applications  ranging  from  semiconductors  to  machine  tools.  A  notable 
derivative  of  the  physical  vapor  deposition  process  is  one  called  ion  plating[2].  In 
the  ion  plating  process,  the  coating  material,  in  atomic  or  molecular  form,  is  to  be 
charged  and  accelerated  in  an  electric  field  so  that  the  charged  particles  impinge  on 
the  substrate  with  a  very  high  kinetic  energy.  The  impact  of  the  charges  particles 
causes  sputtering  of  the  substrate  and  a  great  deal  of  physical  mixing  of  the  coating 
and  substrate  materials.  Usually  what  results  is  a  coating  of  gradual  chemge  in 
material  composition  from  the  substrate  to  the  coating  material.  These  processes  have 
demonstrated  that  composite  materials  with  distinctive  material  properties  but  a 
gradually  varying  interfacial  zone  composition  and,  as  a  result,  varying  interface 
properties  can  be  achieved;  see  for  example  [3].  The  immediate  benefit  of  this 
achievement  is  perhaps  the  relief  of  residual  stresses  in  the  composite.  Because  the 
material  characteristics  at  the  interface  is  no  longer  an  abrupt  change  of  properties 
but  rather  a  smooth  transition  from  one  material  to  another,  thermal  mismatch  and 
residual  stresses  are  reduced  substantially. 

In  another  area  of  the  development  of  composite  materials,  a  concept  that  seeks 
to  design  engineering  composites  to  optimize  material  properties  for  particular 
applications  is  producing  more  new  material  forming  processes.  Prompted  by  the 
"Space  Plane  Project"  that  requires  high  temperature  materials  to  protect  the  space 
vehicle  from  thermal  damage  during  reentry,  one  such  "Functionally  Gradient 
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Materials"  describes  a  successful  laboratory  testing  of  a  TiB2-Cu  system  to  be  used  on 
the  exterior  of  the  spacecraft.  The  TiB2:Cu  ratio  decreases  from  exterior  of  the  surface 
to  the  interior[41. 

In  high  temperature  engine  combustion  chambers  and  gas  turbine  air  foils, 
where  good  heat  resistant  property  near  the  surface,  heat  dissipation  in  the  bulk, 
minimal  residual  stresses  and  overall  mechanical  toughness  are  required,  this  new 
way  of  materials  design  can  succeed  where  traditional  design  methods  fail.  There  are 
several  methods  to  produce  such  "gradient"  effect  in  composite  materials,  such  as 
Chemical  Vapor  Deposition[5],  Centrifugal  Casting[6],  and  Combustion  Sintering 
[7].  Composite  materials  so  produced  possess  mechanical  and  thermal  properties, 
instead  of  changing  abruptly  from  one  to  the  other  as  traditional  composites  do,  that 
vary  continuously  through  the  thickness  of  the  composite[8].  The  resulting 
continuously  varying  material  composition  leads  to  reduced  residual  stresses. 

This  study  is  directed  toward  the  fracture  mechanics  of  an  interface  at  a  two 
material  junction  where  both  materials  are  of  finite  thickness.  The  material  is 
isotropic,  but  near  the  interface  it  is  nonhomogeneous.  This  nonhomogeneity  results 
either  from  intentional  material  mixing  or  through  the  particular  process  used.  To 
study  the  fracture  mechanics  of  the  interface  we  assume  the  existence  of  initial  defects 
in  the  form  of  a  crack,  and  then  try  to  find  certain  fracture  parameters  that  govern 
the  tendency  of  the  crack  to  grow.  In  brittle  fracture  and  fatigue  crack  growth,  the 
most  important  parameter  we  seek  is  the  stress  intensity  factor,  which  combines  the 
effects  of  applied  loads,  crack  geometry  and  material  properties.  To  be  more  specific, 
we  shall  phrase  the  problem  as  follows:  To  study  the  fracture  mechanics  problem  of 
a  crack  at  the  interface  of  a  two  bonded  finite  thickness  materials.  Both  materials 
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will  be  isotropic  one  being  homogeneous  and  the  other  nonhomogeneous  in  the 
thickness  direction.  The  problem  will  be  formulated  as  a  plane  strain  or  generalized 
plane  stress  problem. 

The  analysis  of  this  interface  crack  problem  will  be  performed  within  the 
confines  of  the  linear  theory  of  elasticity,  which  we  solve  as  a  plane  strain  or 
generalized  plane  stress  problem  in  a  straightforward  manner,  under  self  equilibrating 
surface  tractions.  The  crack  problem  is  of  the  mixed  boundary  value  type  and  is 
reduced  to  a  system  of  singular  integral  equations.  Except  under  very  special 
circumstances  where  a  closed  form  solution  can  be  found,  most  often  the  system  of 
SIE’s  are  solved  by  a  numerical  method  to  obtain  the  crack  tip  stress  intensity  factors, 
the  strain  energy  release  rate  and  the  direction  of  probable  crack  extension.  Once  the 
integral  equations  are  solved,  the  stress  field  can  be  readily  computed  if  one  so 
desires. 


1.2  literature  surve 


In  studying  the  fracture  mechanics  of  bonded  materials,  in  recent  past 
considerable  attention  has  been  directed  toward  studying  the  mechanical  behavior  of 
the  interfacial  regions  where  defects  usually  in  the  form  of  voids  or  cracks  often  exist. 
Numerous  publications  have  addressed  the  mechanics  of  interfacial  regions  in  bonded 
materials  (see  for  example  [9]  through  [13]).  It  has  been  shown  that  for  interface 
cracks  the  material  properties,  in  particular  the  ratio  of  Young’s  moduli,  plays  an 
'mportant  role  in  their  fracture  mechanics  behavior.  In  the  simplest  case  of  two 
bonded  materials  where  both  are  homogeneous  and  isotropic,  the  Young’s  moduli  are 
constant,  and  therefore  a  jump  discontinuity  exists  at  the  interface.  With  such  a 
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discontinuity  it  is  well  known  that  the  stress  state  around  the  crack  tip  exhibits  an 
oscillation[9][10][ll].  This  is  physically  inadmissible,  as  it  implies  wrinkling 
and  overlapping  of  the  crack  surfaces  near  the  crack  tips[12].  It  could  be  argued 
that  this  affects  a  very  small  region  near  the  crack  tip  and  consequently  is  not 
important. 

The  crack-contact  model  tries  to  resolve  this  inconsistency  by  assuming  that 
the  crack  surfaces  close  near  the  crack  tips  and  form  a  cusp  and  the  shear  stress  in 
the  contact  region  is  zero[13].  However,  the  general  belief  is  that  the  resolution 
to  this  physical  anomaly  has  to  come  from  a  more  realistic  modeling  of  the  interfacial 
region  of  the  bonded  materials.  Research  has  shown  that  this  stress  oscillation 
disappears  as  long  as  the  material  property  is  continuous  around  the  crack  tip 
[9][14][15][16].  For  the  interface  cracks  there  is  a  model  that  treats  the 
interfacial  region  as  a  third  nonhomogeneous  material  with  steep  property  gradients 
thus  eliminating  the  jump  discontinuity  in  the  material  properties[16].  This  model 
eliminates  the  stress  oscillations  since  the  material  properties  are  continuous 
throughout  the  domain. 

In  light  of  the  discussions  in  the  previous  section  on  composites  that  possess 
varying  material  properties  near  the  interface,  material  nonhomogeneity  is  something 
that  has  to  be  taken  into  account  in  the  study  of  interface  crack  problems  in  a  wide 
variety  of  applications.  Reference  [16]  exemplifies  such  a  case  by  assuming  one 
material  to  be  nonhomogeneous  with  the  material  parameters  varying  through  the 
thickness  in  a  certain  exponential  manner  and  provides  the  interface  crack  solution 
for  two  bonded  half  planes.  In  many  practical  applications  where  the  crack  size  is 
small  as  compared  to  the  thickness  dimension,  often  the  bonding  of  two  materials  can 
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be  modeled  by  two  half  planes.  In  some  other  applications,  however,  this  may  be 
realistic.  Thus,  in  cases  such  as  the  dry-film  lubricant  coatings  discussed  earlier,  thin 
layer  of  refractory  coatings  on  cutting  tools  and  variety  of  components  involving  thin 
films  in  microelectronics  where  the  film  thickness  is  sometimes  less  than  1.0pm  (10"* 
meter)  [3],  it  may  be  necessary  to  investigate  the  interfacial  crack  problem  for  finite 
domains. 

1.3  Overview 

The  main  objectives  of  this  study  are  to  study  the  interface  crack  problem  in 
a  two-layered  solid  with  finite  thickness.  The  medium  is  nonhomogeneous,  and 
particular  attention  is  focused  on  thin  film  cases  where  numerical  computation  is 
difficult.  The  study  admits  any  general  loading  configurations,  covers  a  very  wide 
range  of  material  nonhomogeneity,  and  does  not  place  limits  on  the  material  thickness 
to  crack  length  ratio.  However,  if  the  material  is  strongly  nonhomogeneous  and  at  the 
same  time  the  thickness  to  crack  length  ratio  very  small,  the  effort  in  numerical 
computation  would  increase  in  a  dis-proportionate  manner. 

In  Chapter  2  the  problem  is  formulated  under  considerations  of  plane  elasticity 
by  using  the  two  displacement  components  as  the  primary  unknowns.  Through 
integral  transform  the  problem  is  transformed  to  solving  a  pair  of  coupled  second 
order  ODE’s.  With  the  introduction  of  yet  another  new  set  of  unknowns,  i.e.  the 
derivative  of  the  crack  opening  displacements,  and  applying  the  mixed  boundary 
conditions,  for  the  plane  problem  the  original  formulation  is  reduced  to  a  system  of 
singular  integral  equations  with  the  derivative  of  the  crack  opening  displacements  as 
the  unknown  functions.  The  pair  of  singular  integral  equations  are  of  the  first  kind 
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and  both  have  a  simple  Cauchy  kernel.  As  a  result  the  stress  field  around  the  crack 
tips  exhibits  the  well-known  square-root  singularity. 

Chapter  3  describes  the  method  to  solve  the  pair  of  SEE’s,  namely  the  step  by 
step  procedure  of  transforming  the  SIE’s  to  a  system  of  linear  algebraic  equations,  the 
solutions  of  which  are  the  finite  part  of  the  density  functions.  Of  particular  interest 
is  the  numerical  computation  scheme  to  compute  the  infinite  oscillating  integrals  for 
the  evaluation  of  the  related  Fredholm  kernels.  Special  attention  is  paid  to  the  case 
where  the  material  is  strongly  nonhomogeneous,  and  the  material  thickness  to  crack 
length  ratio  is  small,  and  as  a  consequence  the  convergence  of  the  solution  is  slow. 
An  important  point  that  is  addressed  in  Chapter  3  involves  the  integration  of  a  sign 
(step)  function  and  a  logarithm  function.  Failure  to  pay  special  attention  to  these 
integrals  deters  convergence  of  the  numerical  scheme. 

In  Chapter  4,  a  particular  case  of  very  large  thickness  to  crack  length  ratio  is 
computed  so  as  to  compare  with  known  results.  This  serves  as  a  benchmark  of 
verification,  even  though  the  two  problems  are  not  exactly  the  same.  General  cases 
considered  to  be  of  practical  use  is  computed  under  the  opening  mode  and  shear  mode. 
Also,  the  results  are  discussed  and  interpretation  of  the  trends  of  the  crack  tip 
characteristics  as  a  function  of  the  nonhomogeneity  parameter  and  geometry  changes 
is  examined.  Chapter  5  gives  the  conclusions  and  points  to  direction  of  future 
research. 

Hidden  in  the  straightforwardness  of  the  derivation  in  Chapter  2  is  the 
analytical  formulation  of  the  interface  crack  problem.  The  complexity  of  the 
expressions  of  the  Fredholm  kernels  in  terms  of  the  basic  variables  such  as 
thicknesses,  Poisson’s  ratio  and  the  nonhomogeneity  parameter  can  be  seen  in 
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Appendix  A.  In  it,  all  the  elements  of  the  determinants  that  form  the  numerators  and 
denominators  of  the  expressions  in  Fredholm  kernels  are  listed.  These  expressions 
were  first  derived  using  the  symbolic  manipulation  program  MACSYMA[17]  on  a 
VAX-8300  and  later  verified  by  similar  program  MAPLE[18]  on  a  VAX-8530. 

Appendix  B  describes  the  asymptotic  expansion  of  the  rational  polynomials 
where  the  numerators  and  denominators  are  determinants  the  elements  of  which  are 
depicted  in  Appendix  A.  In  theory,  only  the  leading  terms  need  to  be  extracted  in 
order  to  derive  the  pair  of  SIE’s  with  the  proper  Cauchy  kernel,  which  is  to  be 
expected  for  the  problem  defined.  For  the  sake  of  numerically  solving  the  SIE’s  in  a 
more  efficient  manner,  the  asymptotic  expansion  is  carried  out  to  as  many  terms  using 
the  symbolic  manipulator  MAPLE,  as  the  limit  of  computer  resources  would  allow. 
To  tackle  potential  practical  cases  applicable  to  this  work,  which  usually  involves  thin 
films  deposited  on  substrates  with  film  thicknesses  in  the  micron  range,  the  successful 
asymptotic  analysis  to  as  many  terms  as  possible  proved  to  be  absolutely  essential. 

Appendix  C  gives  some  useful  formulas  to  evaluate  those  infinite  integrals, 
whose  integrands  have  been  known  in  closed  form  through  asymptotic  expansions 
carried  out  in  Appendix  B.  Appendix  D  illustrates,  by  way  of  simple  examples,  several 
schemes  to  numerically  integrate  a  step  function  and  logarithm  function  both  of  which 
have  a  discontinuity  within  the  interval  of  integration.  It  explains  the  technique  used 
in  Chapter  3  to  separate  the  integrals  involving  step  function  and  the  logarithm 
function  from  other  infinite  integrals.  Failing  to  do  so  would  very  nearly  guarantee 
the  breakdown  of  convergence  in  the  numerical  scheme. 

Appendix  E  addresses  two  aspects  that  are  very  important  in  the  numerical 
computation  in  this  work.  One  is  the  limit  of  floating  point  calculations  that  can  be 
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carried  out  for  the  Fredholm  kernels,  which  depends  on  the  range  of  a  real  number 
defined  in  FORTRAN  language  resided  in  the  computer  and  used  to  perform  the 
computation.  The  other  is  how  such  cases  involving  thin  films  on  thick  substrates  are 
to  be  successfully  tackled  within  the  numerical  scheme. 

Appendix  F  shows  the  derivation  of  the  Cauchy  integral  which  is  necessary  to 
study  the  singular  behavior  of  the  solution  around  the  crack-tips  and  to  compute  the 
stress  intensity  factors. 
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Chapter  2 

Formulation  Of  The  Problem 


2.1  Introduction 

The  crack  problem  under  consideration  is  shown  in  Figure  1.  Two  plates  are 
bonded  together  and  it  is  assumed  that  there  is  a  crack  at  the  interface.  It  is  a  plane 
strain,  isotropic,  but  nonhomogeneous,  elasticity  problem.  However,  the  derivation 
in  what  follows  is  also  good  for  a  generalized  plane  stress  problem  should  any  physical 
interpretation  warrant.  Material  1  (occupying  y  <  0)  is  homogeneous,  whereas 
material  2  (occupying  y  >  0)  is  nonhomogeneous  in  the  y-direction.  The  elastic 
properties  of  material  2  may  vary  very  steeply  near  the  interface  but  continuity  of  the 
material  parameters  is  maintained  across  the  interface.  The  nonhomogeneity  of 
material  2  is  assumed  to  be  such  that  its  shear  modulus  varies  across  the  thickness 
as  an  exponential  function.  This  assumption  is  broad  enough  to  accommodate  most 
practical  applications.  The  real  advantage  made  possible  by  this  assumption, 
however,  is  that  it  leads  to  a  system  of  differential  equations  with  constant 
coefficients.  Given  the  complexity  of  the  nonhomogeneous  problem,  this  feat  is  vital 
in  the  formulation  of  our  problem.  As  for  the  Poisson’s  ratio  v  ,  we  shall  consider  it 
constant. 

To  solve  the  linear  elasticity  problem  under  general  loading  conditions  defined 
above,  we  use  the  method  of  superposition  as  shown  in  Figure  2.  First  the  elasticity 
problem  of  the  configuration  as  defined  in  the  absence  of  any  cracks  under  a 
prescribed  loading  is  formulated,  Figure  2(b).  The  crack  problem  as  shown  in 
Figure  2(c),  the  solution  of  which  is  of  the  main  interest  in  this  study,  is  solved  by 
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applying  loadings  at  the  two  crack  surfaces  that  is  equal  in  magnitude,  but  opposite 
in  sign,  to  the  tractions  at  the  same  locations  of  the  elasticity  problem  in  Figure  2(b). 
Since  only  the  crack  problem  (Figure  2(c))  carries  the  characteristic  stress 
singularities  that  we  seek,  the  solution  of  it  alone  will  tell  us  the  fracture  mechanics 
characteristics  of  the  problem.  Only  when  we  seek  to  know  the  stress  field  need  we 
superpose  the  solutions  of  the  two  elasticity  problems. 


2.2  Formulation  of  the  crack  problem 

Referring  to  Figure  1,  let  the  shear  modulus  of  material  1  be  pj,  which  is  a 

constant,  and  that  of  material  2  be  p^  =  Pj  eYy.  The  thicknesses  of  the  two  materials 

are  h}  and  h2,  respectively.  We  use  the  two  components  of  displacement  u  and  v  as 

the  primary  unknowns.  Let  k  =  3  -  4v  for  plane  strain  and  k  »  — _ X.  for 

1  +  v 

generalized  plane  stress,  where  v  is  the  Poisson’s  ratio.  Hooke’s  Law  may  then  be 
expressed  as 


o„  -  JL[(k*1)^L  +  (3*kA, 
K  - 1  ox  dy 


-  JLj(3-k)4£l  +  (Kii)i 

dy 


w  K-r '  dx 


du  dv 
dx  dy 


_  ,  ou  ov , 

■  ft  (—  +  — ). 


(1) 


i  « 


In  the 


1,  2. 

absence  of  body  forces,  the  equations  of  equilibrium  are  expressed  as 
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Figure  1.  A  crack  in  the  interface  of  two  materials  with  finite  thicknesses. 
One  material  is  homogeneous  and  another  non-homogeneous. 
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(X,h2)  "Oj  Oyy  (XjO-)  —Oyy  (XjO  +  ) 
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Figure  3.  Boundary  conditions. 


(2) 


9c„  9o_ 


dx  dy 


0, 


♦  i3t  -  0. 

dx  dy 

Substituting  Eqn.  (1)  into  Eqn.  (2)  for  0  <  y  <  h2  we  find  the  Navier’s  equations  as 


(K+i)^i+(K-i)i?ii+2^L+7(K-i)iJi+7(K-i)iiL  -  o, 

9a:2  9y2  dxdy  dy  dx 


(k  -  l)i!iL  +  (K  + 1)151  +  2_^_  +7(3  -k)5“  +**  ♦  1)  -  0. 

9x2  9y2  dxdy  dx  9y 


(3) 


Equation  (3)  and  similar  equations  obtained  for  material  2  must  be  solved  under  the 
following  boundary  and  continuity  conditions  (Figure  3): 


GyyixJiJ  -  0,  axy(x,h2)  -  0, 

GyyiXj-hJ  -  0,  o^x,-/^)  -  0,  -oo<x<<», 

oJy(x,0+)  -  o^x.O-),  o^(x,0+)  -  o^x.O-),  — °°<x<°o> 

aJ5,(x,0+)  -  o^O-)  -  -pjx),  cxy(x,0+)  «  o^fx.O-)  -  -p2(x),  -a<x<a, 
u(x,0+)  -  u(x,0-),  u(x,0+)  -  v(x,0-),  x£-a,  xka, 


(4) 

(5) 

(6) 

(7) 


where  p/x)  and  p/x)  are  the  normal  and  shear  stresses,  respectively,  at  the  location 
of  the  crack  in  the  elasticity  solution  of  the  corresponding  uncracked  problem 
(Figure  2(a)).  The  notations  0+  and  0*  are  necessary  since  there  is  a  discontinuity  of 
displacement  field  at  the  crack  surface. 

To  solve  Eqn.  (3),  define  the  Fourier  transforms  of  the  two  displacement 
components,  u  and  v,  respectively,  as  follows: 
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U( ay)  •  J  u(xy)  e',az  dx, 

V(a,y)  *  J"  v(xy)  e"“  dx. 

By  definition,  u  and  v  become 

u(xy)  -  JL  f"  U(ay)eia*  da, 
2nJ— 


v(xy)  -  JL  f“  V(a«y)eio*  da. 
2jc*'— 


(8) 


From  Eqns.  (3)  and  (8)  it  can  be  shown  that 

(K-l)^+Y(K-l)4^-(K  +  l)a2t7+2ia.^+Y(K-l)iaF  -  0, 
dy2  dy  dy 


(K  +  l)^!Z+Y(K  +  l)-^  +  (K'1)a2v+2ia— +Y(3-K)iaf7  -  0. 
dy2  dy  dy 


Solving  Eqn.  (10)  we  find 


Way)  -  C1(a)e",y +C2(a)e"*>’  +C3(a)e"*>’ +C4(a)e',*y, 


(10) 


V(ay)  -  D1(a)eBJ'+D2(a)e’,*5,+D3(a)e'v+I>4(a)en-y, 


(11) 


0  <y  <h2, 


where,  nJt  n2,  n3  and  n4  are  solutions  of  the  following  characteristic  equation: 


n4  +  2Yn3-(2a2-,f)n2-2a2Yn  +  ^_l!ia2Y2+a4  -  0. 

1  +  K 


(12) 


Eqn.  (12),  a  polynomial  equation  of  the  fourth  degree,  has  the  following  roots: 
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The  above  derivation  is  for  the  nonhomogeneous  material.  For  material  1,  which  is 
homogeneous  (y  =  0)  or  for  -h2  <  y  <  0,  Eqn.  (10)  becomes 
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(K-l)^!^-(K  +  l)a2L7+2a^X  -  0, 
dy2  dy 


(K  +  l)f^ -(K-l)a2V+2ia—  -  0. 
dy 2  dy 


The  corresponding  characteristic  equation,  which  can  be  obtained  by  letting  y  =  0  in 
Eqn.  (12),  becomes 


(n2-<x2)2  -  0. 


From  Eqns.  (17)  the  solution  of  Eqn.  (16)  may  be  obtained  as 


U(ay)  -  [A1(a)+A2(a)y]e|ab’+[A3(a)+A4(a)y]e*lob', 
V(a,y)  -  [B1(a)+B2(a)y]elo|,  +  [B3(a)+B4(a)y]e‘ia|>, 


-hx  <  y  <  0. 


Similar  to  the  nonhomogeneous  case,  the  functions  A/a)  and  B/a),  ( j  =  1  — -  4)  are 
inter-dependent.  Their  relationships,  obtained  by  substituting  Eqn.  (18)  into  either 
one  of  Eqn.  (16),  are  found  as  follows: 


B/a)  -  -*[-fLA1(a)-JlA2(a)], 
ia  a 


BJa)  -  -i-2—A2( a), 


B/a)  -  i  [_EL As (a) + JSa4  (a) ], 

lot  I  a 


B4(a)  -  t..a,-A4( a). 

I«  | 
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The  original  problem  has  only  two  unknowns,  i.e.  the  two  displacement 
components.  After  Fourier  transform,  we  have  instead  eight  unknown  coefficients  in 
the  solution  to  the  two  ODE’s.  Now  we  need  to  utilize  the  homogeneous  boundary 
conditions,  Eqns.  (4)  through  (7),  to  determine  the  unknown  functions  A/s  and  C/s. 
Before  we  apply  Fourier  transforms  to  the  boundary  conditions,  define  density 
functions  as  follows: 


fox)  -  JL[i;2(;c,0+)-i;1(x,0-)], 

OX 

fox)  -  JLti^O^O+l-UjOe.O-)]. 

dx 


(20) 


It  follows  that  fox)  =  fox)  =  0,  -oo  <  x  <  -a,  a<  x  <  °°,  which  is  a  direct  result  of  the  fact 
that  v2(x,0+)  =  Vj(x,0-)  and  u2(x,0+)  =  u1(x,0-)  outside  of  the  crack.  Let  the  Fourier 
transforms  of  /z  and  f2  be  Fz  and  F2,  respectively,  that  is,  let 


F,(a)  -  fcfox)  e'iax  dx, 
F2(a)  -  iafox )  e'iax  dx. 


(21) 


The  Fourier  transforms  of  Eqns.  (4)  -  (7)  may  then  be  expressed  as  follows: 


afox,h2)  -  0,  -oo<x<oo, 


ayy(x,h2)  0,  -oo<x<», 


ojx,-^)  «  0,  -»<*<«, 


__(a,A2)+iaV(a,)i2)  - 
dy 

(3  -K)iaU(.a,h2)  +(k  +  1) 


0, 

dV 

dy 


(a  ,h2) 


——(a, -/ijl  +  iaVfa, -h2)  -  0, 
dy 


0, 


(22) 

(23) 


(24) 
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(25) 


OvyOe.-Ai)  -  0,  -°o<x<oot 


(3  -K)ia(7(a,-A1)  +(K  +  l)^X(a,-A1) 

dy 


0, 


ow(x,0+)  -  o^Oc.O-),  -oo<x<oo, 


(3  -K)ia£7(a,0+)  +  (K  +  l)^X(a,0+) 

dy 


(3-K)iat/(a,0-)+(K  +  l)^X(a,0-), 

dy 


(26) 


ciy(x,0+)  -  o^^O-),  -oo<x<«,  -» 

i^(a,0+)  +iaV\a,0+)  -  —  (a,0-)  +  iaV(a,0-). 

dy  dy 


(27) 


Note  that  the  Fourier  transforms  of  the  density  functions  are  identically  zero  outside 
of  the  crack  and  are  the  new  unknowns  inside  of  the  crack.  We  can  now  express  the 
transforms  of  the  mixed  and  final  boundary  condition  as 

ia[V(a,0+)  -V(a,0-)]  -  Fj(a), 

(28) 

ia[f/(a,0+)  -  £7(a,0-)]  -  F2(a). 


We  now  substitute  the  appropriate  expressions  for  IT s  and  V’s  from  Eqns.  (11)  and 
(18)  into  Eqns.  (22)  to  (27).  By  taking  advantage  of  the  inter-dependence  between  A/s 
and  BjS,j  =  1  4,  and  C/s  and  D/s,  j  =  1  — •  4,  we  obtain  a  system  of  eight  linear 

equations  for  the  unknown  functions  A/s  and  C/s,  (J =  1  ••••  4)  in  terms  of  the  Fourier 
transforms  of  the  density  functions  as  follows 
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0 

0 

0 

0 

ai5 

ai6 

ai7 

ai8 

Ajfa) 

0 

0 

0 

0 

0 

a2i 

°26 

a27 

a2B 

A2(a) 

0 

a31 

fl32 

a33 

°34 

0 

0 

0 

0 

A3(a) 

0 

a*l 

a42 

a43 

a44 

0 

0 

0 

0 

A4(a) 

0 

fl51 

aS2 

«53 

ab4 

a5S 

aS6 

a57 

a&8 

Cj(a) 

0 

a61 

°62 

a63 

flfi5 

°6€ 

°67 

aS8 

C2(a) 

0 

< hi 

a72 

a73 

Ch4 

a75 

a76 

a77 

®78 

C3(  a) 

Fj(a) 

_  a81 

0 

a83 

0 

aBB 

a8fi 

°B7 

a8S 

_C4(a)  _ 

F2(a) 

(29) 


where  aLJ ,  i,  j  =  1  ••••  8,  the  exact  expressions  of  which  are  given  in  Appendix  A,  are 
functions  of  a,  y,  k,  h2,  and  h2.  We  can  solve  Eqn.  (29)  for  Ait  Ct,  ( i  =  1  •••■  4)  in  terms 
of  the  unknowns  Fp  F2  as 


(30) 


Ci.4  -  E  —Fj>  im  5”‘8> 
j-i  D 


(31) 


where  D  is  the  determinant  of  the  8  by  8  coefficient  matrix,  and  DtJ  are  the 
appropriate  cofactors  of  D.  is  defined  as  the  7  by  7  determinant  computed  from 
deleting  the  i-th  column  and  the  j+6- th  row  of  the  8  by  8  matrix,  then  multiplied  by 
the  factor 


Up  to  this  point,  we  have  used  all  of  the  boundary  conditions  except  the 
traction  boundary  conditions  on  the  crack  surfaces  (Eqn.  (7)),  namely, 
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o^Ot.O-)  -  -PjOt),  c^Oc.O-)  -  -p2(x),  -a<x<a, 


(32) 


or 

<^(*,0+)  -  -pi(jc),  axy(x,0+)  -  -p2(x),  -a<x<a.  (33) 

Either  of  (32)  or  (33)  can  be  used;  we  chose,  however,  (32)  for  its  simpler  expressions 
because  material  1  is  homogeneous. 

We  now  have  two  equations  (32)  in  which  we  substitute  stresses  expressed  in 
terms  of  Fourier  transforms  of  the  displacement  components  by  using  the  Hooke’s 
Law.  These  transformed  displacement  components  are  given  by  (11)  and  (18)  in  terms 
of  eight  undetermined  coefficients  or,  indirectly  in  terms  of  Fourier  transforms  of  the 
two  density  functions  introduced  by  (21).  Therefore,  the  two  equations  (32)  are 
sufficient  to  determine  the  unknowns  Fj  and  F2. 


2.3  Derivation  of  the  system  of  singular  integral  equations 

Expressing  Eqn.  (32)  in  terms  of  Fourier  transforms  of  the  displacement 
components,  we  obtain 


■p1(x)  -  - ^ - [(3  -k)  f- U(<x,0-)iaeiaz  da  +  (x  + 1)  f“ ^X(a,0-)eiO1  da], 

2x(k-1)  dv 


dy 


(34) 


-p2( X)  «  -E!L[(3  -  k)  f~i_(a,0-)eiax  da  +  (K  +  l)f" V(a,0-)eio*  da], 
2n  dy 


By  substituting  from  Eqn.  (18)  into  Eqn.  (34)  it  may  be  shown  that 
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(35) 


Jj(jc)  «  lim  r"-m({2Aj  +  [2y-i(ic  +  l)_J_]A2}  eiaj!'|a!-v  + 
v_  o-  J—  a 


{2A3  +  [2y  +  i  (k  +  1)_L_]AJ  e‘“'la|y)  da, 

I  a  | 

-lEpj(x)  -  lim  f"{[2|a|A1+(2|a|y-K  +  l)Aa]  eiai*|a|*  + 

Pi 

[-2  |a  |A3-(2  |a  \y  +  k-1)A4]  eia*~laly)da. 

Now  substituting  (30)  into  (35),  expressing Fj,  and  F2  in  their  integral  iorrn  (Eqn.  (21)) 
and  then  changing  the  order  of  integration,  we  may  write  th.  integral  equations  with 
the  density  functions  as  unknowns  as  follows: 


-  -  fa  Y,KM)mdt  i  -  1,  2,  -a<x<a, 

2  ^  j. i 

where 

KM)  -  lim  JUlf"  -ia({2^i  +  [2y-(K  +  l)_L]£!i}  + 

^  y-  o-  4  D  jai  D 


{ 2^11  +  [2y  +  (k  +  l)_l_]^Si}  da,  j  -  1,  2, 

D  jai  D 


K2Ax,t)  -  lim  JSlif" ([2  | a  ii-  (2  |a  |y-K  +  l)]^L]  e*'**-”*^- 
0-  4  D  D 


[2|o|^L  +  (2|a|y  +  K-l)]i]  e^-o-i-b)  da,  j  -  1,2. 


It  should  be  pointed  out  that  the  integrands  in  Eqns.  (37)  and  (38)  axe  complex 
expressions  involving  determinants  of  several  7  by  7  and  8  by  8  matrices  the  elements 
of  which  are  also  long  expressions  of  complex  variables,  as  detailed  in  Appendix  A. 
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But  the  boundedness  of  these  integrands  at  a  =  0,  can  be  seen  by  noting  that  the 
integrands  vanish  at  cx  =  0.  We  also  note  that  the  integrands  in  Eqns.  (37)  and  (38) 
are  continuous  functions  of  a.  It  then  becomes  clear  that  any  singularities  the  kernels 
KtJ  must  come  from  the  asymptotic  behavior  of  the  integrands  as  |a|  approaches 
infinity. 

The  details  of  asymptotic  expansions  of  the  kernels  are  given  in  Appendix  B. 
In  Eqns.  (37)  and  (38),  the  behavior  of  Rl  as  a  goes  to  «  and  respectively  are 
as  follows, 


BlUa)  - K.-.L_  *  OfJL), 

D  2a(x  +  l)  a2 


£ii(o)  -  ~L  +  O(-l), 

D  2a  a2 


ElLia)  -  _L  +  0(1), 

D  k  +  1  a 


i_S(a)  -  _ L_  ♦  0(1), 
D  k  +  1  a 


lli(P)  -  .  Jili—  +  0( JL). 

D  2p(K  +  l)  p2 


-JL  +  o(JL), 

2P  p2 


£l(P) 

o 


J_  +  0(1), 

K+1  p 


lf£(P)  -  ~_L  +  0(1), 

D  K  +  1  P 


JJ31 

where  P  =  -a  in  Eqn  (40).  We  note  that,  as  shown  in  Appendix  B,  terms  like 
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-~^32  D*1  and  Eli  are  analytic  at  lot  I  =  »  and  they  are  of  much  lower  order 

D  D  DDDD  D 

as  compared  to  .  -11,  .  _iL,and  and  therefore  no  asymptotic  expansion 

needs  to  be  done. 

Let  us  now  try  to  separate  the  leading  terms  from  Eqns.  (37)  and  (38). 

D 

Substituting  the  first  terms  in  Eqns.  (39)  and  (40)  for  _Ji  ,  i,  j  =  1,  2,  separating 
the  infinite  integral  from  to  »  into  two  parts  at  0  and  making  a  change  of  variable 
for  the  part  from  -»  to  0  by  letting  P  =  -a,  after  some  simplifications  we  obtain  the 
leading  term  of  Kn(x,t)  as  follows: 


lim  f  ( 1  +ay)  eaysina(£  -x)da. 
Jo 


(41) 


0 


It  is  easily  shown  that 


lim  f  eay  sina« -x)da  «  - i - 

y—o-  y2+(<-jc)2 

(42) 

1 

■t  _ 

t-x' 


From  Ref.  [19],  we  obtain 

lim  f  ay  ea,sina(£ -x)da  -  0.  (43) 

,-.0-  Jo 

Thus,  in  (37)  the  Cauchy  kernel  has  been  separated.  The  remainder  of  the  integrand 

after  the  leading  term  of  Eli  and  Eli  has  been  separated  is  headed  by  a  term  no 

D  D 

higher  than  —  .  It  is  analytic  everywhere  and  we  can  eliminate  the  limiting 
a 

process  by  substituting  y  =  0  directly  under  the  integral.  This  gives  the  Fredholm 
kernel  ku(x,t)  as  follows 
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(44) 


kn{x,t)  -  JSlIf"  -t  a(  { 2 (BlL  -  — IL-1 — )  - (k  +  l)_J_(^2i  -  _J_)}  e‘a(*- 
4  D  2a(K  +  l)  |a|  D  k  +  1 


[iElL+iK  +  D—lElL)  e *•*-**)  da. 
D  a!  D 


Likewise,  the  leading  term  of  K2/x,t)  is  found  to  exhibit  the  same  singular 
characteristics  as  Eqn.  (42).  The  other  Fredholm  kernel  k 2/x,t)  is  obtained  from  Eqn. 
(38)  in  the  same  manner. 

By  a  similar  process,  the  leading  term  for  K12(x,t)  and  K21(x,t)  has  the  following 
form:  (see  [19]) 


lim  (  ay  ea>cosa(x -t)da  -  0. 

y-0- 


(45) 


Therefore,  we  do  not  need  to  separate  the  leading  term  from  KI2(x,t)  and  K21(x,t)  and 
the  Fredholm  kernels  k];/x,t)  and  k2I(x,t)  are  formulated  in  an  analogous  manner  as 
kn(x,t)  done  previously. 

From  the  discussions  above,  we  find  that  the  singular  behavior  of  the  kernels 
come  solely  from  Ku  and  K22  and  we  can  write  Eqn.  (36)  as  a  system  of  singular 
integral  equation  as  follows 


-"l-fiW  ■  f"  £?—-dt+  (aku(x,t)fj{t)dt  +  fak12(x,t)f2(t)  dt, 
2px  t-x  J-a  J-o 

_n(K  +  l)pjx)  .  f  fjt£idt  +  (ak2l(x,t)f1(t)dt  +  (a k22(x,t)f2(t)  dt, 
2u,  J-a  t-x  •>-»  J-o 


(46) 


-a<x<a. 

kjj,  k12>  k21,  k22  are  derived  from  the  kernels  in  Eqns.  (37)  and  (38),  as  described 
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previously.  They  may  be  expressed  in  the  following  form 


ku{x,t)  ■  J^D^a)  sina(£ -x)  da, 


(47) 


k12(x,t)  ■  J*o  Di2(a)  cosa(t  — jc)  da. 


(48) 


k2l(x,t)  »  f  Z)21(a)  cosa(t  -x)  da, 

J  0 


(49) 


k22(x,t)  «  J~Z)22(a)  sinatt  -x)  da, 


(50) 


where, 


Du(a)  -  -J^[2a(£^^--JilL_)+(K  +  l)&-^li-_i_)], 

2  D  D  2a(K  +  l)  D  D  k  +  1 


(31) 


Z?;2(a)  -  -^llf[2a(^i  +  ^)-(K  +  l)(^i-^i)], 


(52) 


D2»  -  Jili[2a(^i-^i)-(K-l)(^i+^i)], 


(53) 


DU  a)  -  -JElii[2a(£ii-fl«-  J_) -(k-1)(Z^ -_L)]. 

2  D  D  2a  D  D  k  +  1 


(54) 


We  now  have  derived  the  system  of  singular  integral  equations  (Eqn.  (46))  with 
a  Cauchy  kernel.  The  SIE’s  will  be  solved  for  the  density  functions  f}  and  f2  which  are 
the  unknowns.  There  are  two  additional  conditions,  namely, 


f*  ftt )  dt  -  0,  i  -  1,  2, 


(55) 


which  are  the  so  called  single-valuedness  conditions.  Physically,  they  mean  that  the 
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crack  opens  at  one  end  and  closes  at  the  other.  We  will  show  that  with  Eqns.  (46)  and 
(55),  we  can  solve  the  SIE’s  numerically  by  transforming  them  to  a  system  of  linear 
algebraic  equations. 
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Chapter  3 

Solution  Of  Singular  Integral  Equations 
And  Numerical  Computations 


3.1  Introduction 

In  Chapter  2,  we  have  derived  the  system  of  singular  integral  equations  for  the 
interface  crack  perturbation  problem  for  finite  thickness  materials.  These  SIE’s  each 
has  a  Cauchy  kernel  and  two  Fredholm  kernels.  The  two  unknowns  to  be  solved  are 
the  density  functions  which  are  the  derivatives  of  the  difference  of  the  two 
displacement  components.  The  solutions  of  SIE’s  are  generally  obtained  either 
through  function  theoretical  technique  as  given  by  Muskhelishvili  in  [20],  or 
through  numerical  methods  [21],  [22].  For  the  problem  at  hand  the  SIE’s 
which  have  a  simple  Cauchy  kernel  and  lengthy  Fredholm  kernels,  it  is  most 
convenient  to  use  a  numerical  method  to  obtain  the  solution.  In  this  work  the  method 
of  collocation  with  the  unknown  functions  represented  by  Chebyshev  polynomials 
described  in  [21]  is  used. 

3.2  Solution  of  the  integral  equations 

To  facilitate  solving  the  integral  equations,  let  us  make  a  change  of  variables 
as  follows: 

t  -  sa,  x  -  ra,  =>  dt  -  ads,  dx  -  adr  (56) 

The  singular  integral  equations  (Eqn.  (46))  then  become 
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3.2.1  The  infinite  integrals 

To  evaluate  the  Fredholm  kernels  which  contain  integrals  with  infinite  upper 
limit  (Eqns.  (58)  through  (61)),  we  shall  treat  the  integrands  with  sine  and  cosine 
terms  separately.  The  integrals  that  contain  sine  functions  as  in  Eqn.  (58)  may  be 
expressed  as 
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(62) 


kn(r,s)  -  J^A(£>1‘1(a)-—L)sinaa(s-r)  da  + 

f  (01*1(a)-_JL)sinaa(s -r)  da  +  2.  f  sinaa(s  r^a) 
ja  4a  4  Jo  a 

where  the  integrand  of  the  first  integral  is  bounded  everywhere  and  is  integrated 
numerically.  The  third  integral  has  a  closed  form  expression  given  by 

f~-L  sinaa(s-r)  da  m  sign[a(s-r)]l 1,  (63) 

Jo  a  2 

where  the  sign  function  is  defined  as  a  function  that  gives  only  the  sign  of  its 
argument  with  the  numerical  value  1.0.  The  second  integral  can  be  written  as 

I  (£>,’,(  a)-_I-)sinaa(s  -r)  da  -  |(d,',(a)--jL)sinaa(s-r)  da  + 

Ja  11  4a  Ja  11  4a 


i !(a)]sina a  (s-r)  da. 


Let  0( a")  denote  polynomials  of  degree  n  in  a.  Note  that  du'( a)  is  a  polynomial  of 

finite  degree  in  a  shown  in  Appendix  B.  It  represents  the  leading  asymptotic  terms 

of  Du'(a)  and  the  last  term  in  the  polynomial  is  a  term  like  (Jl)11  0(0)  (see  Eqn. 

a 

(B.59)).  Therefore  the  leading  term  in  the  integrand  in  [D}1'(oQ  *  du’(a)]  is  a  term  like 

(JL)12  0(0)  .  Since  y  is  the  parameter  that  determines  the  degree  of 
a 

nonhomogeneity  in  material  2  and  is  a  constant,  we  can  choose  the  value  "A”,  which 

separates  the  infinite  integral  into  two  integrals,  such  that  (J-)12  is  small  to  any 

A 

order  we  desire  so  that 


fcu>h(a)-d  ^(a^sinaais-r)  da 
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is  negligible. 

A  few  words  must  be  said  here  regarding  the  selection  of  the  value  "A"  such 

that  Eqn.  (65)  can  be  neglected.  While  it  is  true  that  we  can  choose  "A"  so  that 

(-I)12  is  small  to  any  degree  we  want,  for  a  particular  nonhomogeneity  constant, 
A 

there  is  a  point  where  any  further  increase  in  "A"  will  only  serve  to  tax  the  numerical 

effort,  whereas  the  final  solution  has  become  stabilized.  The  increased  numerical 

effort  that  comes  as  a  result  of  the  increase  in  "A”  can  be  seen  from  the  first  integral 

in  Eqn.  (62).  Since  that  integral  is  to  be  evaluated  by  Gauss’  formula,  a  larger  upper 

limit  of  integration  means  more  computing  effort. 

When  a  proper  "A"  has  been  chosen  in  Eqn.  (62),  Eqn.  (64)  can  be  evaluated  in 

closed  form  since  (d^a)  -  —L~)  is  a  polynomial  of  finite  degree  in  a.  The  closed 

4a 

form  expressions  for  the  integral  are  derived  in  Appendix  C.  Note  that 


li 


dija)  -  _L  -  £  On  JL, 

4a  ft  Uj  a> 


(66) 


where  cnJ*,  j  =  2  •••<  11  are  shown  in  Eqn.  (B.59)  in  Appendix  B.  Eqn.  (62)  may  now 
be  written  as  follows 


kn(r,s)  -  \  (D[Aa)  -  _L)sinaa(s  -r)  da 
Jo  4a 


(67) 


ii 


£  .  f-  sinaKs-r)  da  ,  1  sjgna(s_r)5 

ft  Ja  a>  4  2 


Similarly,  k22(r,s)  may  be  expressed  in  the  following  manner: 
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k22(r,s )  -  \  (D22(ct)  -  _I_)8ina a(s-r)  da  + 

Jo  4a 


(68) 


f  .  r-  smaa(s-r)  da  f  y  signa(s_r)*. 
fZ  J  a/  4  2 

We  see  that  k11(r,s)  and  k22(r,s)  are  both  bounded  functions,  even  if  (s  ■  r)  goes 

to  zero.  In  Eqns.  (67)  and  (68),  only  the  first  term  is  to  be  integrated  numerically;  the 

second  term,  which  is  expressed  in  summation  of  integrals,  can  be  expressed  in  closed 

form.  The  closed  form  expressions  for  these  integrals  are  shown  in  Appendix  C.  The 

last  term,  however,  is  a  step  function.  As  (s  -  r)  goes  from  positive  to  negative,  the 

value  of  the  step  function  changes  from  HL  to  -JH.  While  the  reason  for 

8  8 

separation  of  this  term  from  the  Fredholm  kernels  kn(r,s)  and  k2/r,s)  may  not  be 
immediately  apparent,  the  fact  remains  that  there  is  this  term  hidden  in  the  kernel. 
We  shall  illustrate  by  way  of  examples  in  Appendix  D  that  failure  to  take  special 
precautions  when  dealing  with  functions  that  contain  a  step  function  in  numerical 
evaluation  of  integrals  will  lead  to  erroneous  results,  thereby  inhibiting  convergence 
of  the  final  solution.  Furthermore,  the  definite  integral  (see  Eqn.  (57)) 

.15  signals -r)  fx(s)  ds,  i  •  1,  2.  (69) 

of  the  step  function  has  a  trivial  closed  form  expression  as  will  be  shown  later  in  the 
chapter.  Had  the  step  function  not  been  separated,  we  will  find  that  numerical 
convergence  would  be  greatly  hindered  simply  because  of  the  inability  of  numerical 
quadrature  to  properly  approximate  the  step  function. 

Having  described  the  way  to  evaluate  kjj(r,s)  and  k2/r,s),  we  now  do  the  same 
for  kj/r,s)  and  k21(r,s).  The  kernel  k2/r,s)  may  be  expressed  as 
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k12 (r,s)  ■  J^Ar>i2(a)cosaa(s -r)  da  + 

(70) 

f"(D^a)-JL)cosoa(«-r)  da  *  1  f'C08aa(s-r)da. 

Ja  4a  4.Ja  cl 

Again,  the  integrand  of  the  first  integral  in  Eqn.  (70)  is  bounded  everywhere  within 
the  limits  of  integration.  The  third  integral  has  a  closed  form  expression  as  follows 
(see  also  Appendix  C) 

f-  do  •  -CUAals-r » 

Ja  a 

(71) 

1  i  A  /  .  j  /’lAo(*-r)|  C0Sa  ~  1  . 

-  -y0  -  log  jAa(s  -r)  |  -  — - da, 

Jo  a 

where  y„  =  0.57721566490,  is  the  Euler’s  constant.  By  following  the  same  argument 
as  used  in  evaluating  the  second  integral  in  Eqn.  (62),  we  can  choose  a  sufficiently 
large  "A”  so  that  (DI2(aJ  -  d12(a))  is  small  to  any  order  we  desire,  which  reduces  the 
second  integral  in  Eqn.  (70)  to 

(  (D1*2(a)-_L)cosaa(s-r)  das  |  (di2(a)-_L)cosaa(s-r)  da 

Ja  4a  Ja  4a 

(72) 

ii  *  , 

»  52  cizj  —  cosaa  (s-r)  da. 

j. 2  jA  a/ 

Each  term  of  the  summation  in  Eqn.  (72)  is  evaluated  in  closed  form,  similar  to  Eqn. 
(64).  Therefore,  the  Fredholm  kernel  k, /r,s)  may  be  written  as 
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(73) 


ku(r,s)  •  J*  D^oOcosaoKs -r)  da  + 

£c-  f-  ™?«a(s-j2  da  -  IatA«<«-r», 
j-2  jA  o'  4 

where  Ci  is  defined  in  Appendix  C.  The  reason  that  Eqn.  (71)  is  separated  from 
k12(r,s)  is  similar  to  why  the  step  function  is  separated  from  the  kernel  kn(r,s).  It  is 
because  Ci(A  a  ( s-r )),  which  contains  a  logarithmic  term,  could  not  be  properly 
approximated  by  numerical  quadrature  when  integrated.  This  fact  will  be  illustrated 
in  Appendix  D  by  way  of  examples  to  show  the  significant  error  encountered  when 
using  Gauss’  formula  to  numeric  ally  integrate  a  logarithmic  function.  Therefore  the 
logarithm  term  is  separated  to  be  integrated  in  closed  form  as  will  be  seen  in  section 
3.2.3. 

Similarly,  k21(r,s)  may  be  written  as 
k21(r,s)  *  •D2*1(a)cosaa(s-r)  da  - 

(74) 

£  .  r-  cosa  a  (s-r)  ^  t  rCilAa(s-r». 
j.  2  Ja  a7  4 

It  is  worth  noting  that  k12(r,s)  and  k21(r,s)  differ  only  by  a  sign,  as  is  the  case  for 
problems  that  possesses  symmetry  with  respect  to  the  y-axis. 

3.2.2  Nature  of  the  stress  singularity 

It  is  well  known  in  fracture  mechanics  that  the  stress  field  near  the  crack  tip 
possesses  a  characteristics  that  is  proportional  to  rp,  where  r  is  a  small  distance 
measured  from  the  crack  tip  where  we  evaluate  the  stresses,  p  is  called  the  power  of 
singularity,  a  measure  of  the  unboundedness  of  the  stress  field,  whose  value  should 
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lie  between  0  and  1.  If  p  is  larger  than  1,  the  displacement  field  and  the  energy  are 
both  unbounded  at  the  crack  tip,  which  is  unacceptable.  On  the  other  hand,  if  p  is 
smaller  than  0,  r'p  maintains  a  positive  power  and  the  stresses  vanish  as  r  approaches 
0  and  there  is  no  stress  singularity.  In  crack  problems,  we  accept  bounded 
displacements.  The  stresses,  which  are  of  the  order  of  one  less  than  the 
displacements,  may  be  unbounded,  but  must  be  integrable. 

From  the  derivation  of  the  singular  integral  equation  in  Chapter  2  and  the 
discussion  earlier  in  this  chapter,  we  know  that  the  terms  in  Eqn.  (57)  that  involve 
Fredholm  kernels  are  bounded  and  analytic  everywhere.  It  is  only  the  terms  with 
Cauchy  kernel  that  have  singular  behavior.  Furthermore,  from  the  discussions  in 
Appendix  F  concerning  the  characteristics  of  the  density  functions  f/s),  f2(s),  we  note 
that  they  are  of  the  form: 


fSs) 


Ft(s) 


F^sloKs)  ,  i  -  1,  2. 


(75) 


(s  +  l)w(l-s)w 

where  F/s),  i  =  1,  2  satisfy  the  Holder  condition  on  the  closed  interval  of  [-1,  I]  and 
Ffl)  *  0,  F/l)  *0,i  =  1,2. 


3.2.3  Converting  to  linear  system 

Having  described  the  way  these  Fredholm  kernels  are  evaluated,  and  the 
nature  of  the  unknown  functions,  we  are  now  ready  to  solve  the  singular  integral 
equations  by  converting  them  into  a  linear  system  of  algebraic  equations  to  obtain  the 
unknown  functions  at  discrete  collocation  points. 

Let  us  express  our  unknown  functions  in  the  singular  integral  equation  in 
terms  of  Chebyshev  polynomial.  T/s)  and  U/s)  are  the  Chebyshev  polynomial  of  the 
first  kind  and  second  kind  respectively,  and  they  are  defined  as  follows 
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(76) 


Tj(s)  -  cos/0, 

Uj(s>  -  sin(-'>1)e 
sin0 


s  ■  COS0. 


F;(s)  -  Tj(s),  i  -  1,2.  (77) 

j'O 

are  now  the  new  unknowns  to  be  solved  in  the  linear  system.  Substitute  Eqn.  (77) 
into  Eqn.  (57),  keeping  in  mind  the  following  identity, 


2  p  Tft)  dt 
71  ^  ( t-x )  sPT?- 


j  -  o, 


j  >  0, 


Eqn.  (57)  becomes 


-  52  Uj-i(r)  +  —  [  Pfcn(r,s)  fx(s)  a  ds 

2pi  j.  i  rt  (79) 

+  /^(s)  a  <*s],  -l<r<l. 

--lil/j2(r)  -  52B2j  Uj.Jr)*—  [  pA2i(r,s)  fx(s)  a  ds 
2»i  M  (80) 

+  J^fe22(r,5)  f2(s)  a  ds],  -l<r<l. 

Let  us  evaluate  the  integrals  in  Eqns.  (79)  and  (80)  separately.  Substituting  Eqn.  (67) 
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into  the  first  integral  in  Eqn.  (79),  we  obtain 


1  ku{r,s)  f^s)  a  ds  •  & v  a  P  (r,s)  J  S—  ds  + 

1  j-o  -1  /1-s2 

JP  signa(s-r)  f^s)  ds. 


(81) 


fc,*i(r,s)  «  |  (DiAa)  -  — L)  sina a(s-r)  da  + 
Jo  4a 


sinaa  (s-r) 
a* 


(82) 


Recall  in  Chapter  2  the  definition  of  the  density  functions  (See  Eqn.  (20))  and  the 
change  of  variable  we  have  made  at  the  beginning  of  the  chapter 


/i(s)  -  JL  [o2(x,0+)  -  ^(jc.O-)]  a.  (83) 

dx 

The  integral  of  the  second  term  in  Eqn.  (81)  thus  becomes 


£  signa(s-r)/ i(s)  ds  •  J°  sign(t-x)  fjit )  dt 

•  ~jX  /i«)  dt  +  fx{t)  dt 

«  -tu2(i,0+)-i;1(i,0-)]j  +  [u2(*,0+)-i;1(f,0-)]j 

I  -o  I. 

-  -2  [u2(jc,0+)  -  UjU.O-)]. 


(84) 


Because  the  crack  opening  displacement  is  zero  at  the  crack  tips,  Eqn.  (84)  turns  out 
to  be  exactly  twice  the  negative  vertical  crack  opening  displacement  at  point  "x".  This 
crack  opening  displacement  may  be  expressed  as  follows,  after  appropriate  change  of 
variables  x  -  r  a 
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(85) 


v2(x,0+)  -  Vjix, 0-)  -  jT  fyis)  ds 
-  ±BV  £ 

\/l  -s2 


ds 


Again,  make  a  simple  change  of  variable  in  the  integral  in  Eqn.  (85)  by  letting  s  =  cos 
0,  then  ds  =  - sin  0  *29,  and  we  can  easily  show  that  it  becomes 


fr  .  ds  •  -  f°°* V  cos  /  0  *20 

•M  Ft  ~T~  J  * 


n1 

_  sin;  (cos 'V) 

j 

-  -1  Uj.jir)  y/l-r2  . 

j 

Consider  the  Gauss-Chebyshev  integration  formula  [23] 


1  P  .ffilff-is  1  £  git„),  Tn(tk)  -  0, 
*  Vl-*2  n  * -1 


(86) 


(87) 


(„  -  -  X).  k  .  1  •••  n. 

*  2n 


Note  that  to  use  this  formula,  we  need  to  know  the  function  values  of  g( s)  at  points 
tv  k  =  1  ••••  n,  all  zeroes  of  Chebyshev  polynomial  of  the  first  kind  of  degree  n. 
Applying  this  formula  to  the  first  integral  in  Eqn.  (81),  we  obtain 

P  k'n(r,s)  -SfL  ds  -  L  1tk'u(r,ti)  T/tJ.  (88) 

J‘1  n  *-i 

At  any  point  ”r",  and  a  particular  choice  of  zero  of  the  Chebyshev  polynomial  of  the 
first  kind,  the  first  integral  of  kjj*(r,tk  )  (see  Eqn.  (82))  is  evaluated  by  straight¬ 
forward  Gauss’  formula  of  arbitrary  interval  as  follows 
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W.  -  _JL_  [P^)]2.  (90) 

I"*,2 

In  Eqn.  (89),  N  is  the  number  of  sub-intervals  and  it  is  equal  to  the  integer  part  of 

Aid,  q  is  the  number  of  Gauss  points  and  d  is  the  sub-interval  size  for  Gauss*  formula 

evaluation.  As  discussed  earlier  in  Section  3.2.1,  the  choice  of  "A"  is  dependent  on  the 

parameter  y,  the  degree  of  nonhomogeneity  in  material  2,  and  how  small  we  wish  the 

term  ( -X)12  to  be.  This  limit  of  integration  is  therefore  not  a  constant.  The  purpose 
A 

of  separating  the  integral  in  Eqn.  (89)  into  several  sub-intervals  is  to  accommodate  the 
changing  "A"  value  so  that  a  more  uniform  evaluation  using  the  Gauss’  formula  is 
achieved.  Choosing  the  number  of  Gauss  points  is  somewhat  arbitrary.  When  the 
Gauss  points  are  too  few,  we  might  not  have  a  good  representation  of  the  integrand 
hence  an  inaccurate  evaluation  of  the  integral  results.  On  the  other  hand,  it  takes 
only  a  few  trials  to  find  that,  above  a  certain  number,  increasing  the  Gauss  points  will 
only  increase  the  computation  effort  with  no  discernible  increase  in  the  result  of  the 
integral.  For  the  numerical  computation  in  this  work,  a  sub-interval  size  d  -  20  and 
the  number 
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of  Gauss  points  q  =  20  are  used. 


To  complete  the  evaluation  of  Eqn.  (81),  we  need  to  calculate  the  term  that 
involves  integrals  with  an  infinite  upper  limit,  which  we  have  said  in  3.2.1  to  have 
closed  form  expressions.  As  mentioned  earlier,  the  expressions  for  the  integrals  are 
derived  in  Appendix  C,  and  the  coefficients  cn'  are  obtained  and  listed  in  Appendix 
B.  We  write  the  integrals  with  infinite  upper  limit  in  two  parts 


A  .  r-  sinaatt*  -r) 


E  £ 

o-2  A 


da  - 


a" 


(91) 


A  .  f-  sin ao(t-r)  ,  A  .  f-  sinaaft.  -r)  , 

£  c1UB  I  - -1 - da  +  £  Cu;»-i  \  - ~ - da 

n-i  JA  a2n  „.2  JA  a2"  1 


Each  part  may  in  turn  be  written  as 


5 


n  •! 


/; 


sina  {tk  - r ) 


a 


2  n 


da 


E  cY» 

n  ■  1 


[cosA(f*  -r)>  - : - 

U  (2n  -  DLA2"'2-' 


(92) 


sin A(tk-r)Y£ 
j- 1 


(-iy*1(«*-r)2^-1)(2n-2j)! 
(2n  -  DIA2"'2-'*1 


(-1)" 


(tk  -r)2""1 
(2n-l)! 


Ci(A(<t-r))], 
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'11,2" 


"•2 


-1  l 


■-  sinafr*  -r) 


a 


2"  -1 


da  - 


E  c,,.,.,  [cosA«,-r)E - apgggg, - 


n  - 1 

sinA(*4  -r)J2 


(-lV*1^*  -r)2(7_1)(2n  -2/-1)! 
(2n-2)lA2n-2J 


+ 


(93) 


(-1)“* 


l(tk-r)2n'2 
(2n  -2)! 


si(A(tk  -r))]. 


In  summary,  the  numerical  evaluation  of  Eqn.  (81)  requires  terms  from  Eqn. 
(86),  Eqn.  (89)  and  Eqns.  (91)  through  (93),  we  may  write  it  as  follows 

£  kn{r,s)  fj{s)a  ds  -  Tta  52  Bij  l—  E  tE 

j- o  n  *. i  j.i 


wi  (D'n(y,)-—)sinyia{ti-r))  + 
4 


u 

E  c. 

i*2 


•-  sin  a  a(tk-r) 


a‘ 


da  ]  Tj(th)  ♦  T 


0)-i<r> 


vT 


(94) 


where  the  only  integral  in  the  Eqn.  (94)  is  expressed  in  terms  of  summations  in  Eqns. 
(91)  through  (93)  and  all  appropriate  variables  are  defined  earlier  in  this  section. 

Following  the  same  procedure,  the  integral  involving  k2/r,s)  in  Eqn.  (80)  may 
be  expressed  in  terms  of  summations  similar  to  that  involving  ku(r,s)  as  follows: 
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n  N 


P  A22(r,s)  /2(s)a  ds  -  na  Y  B2j  {2  Y  [£) 

1  .-A  W  t  _  1  I  _  1 


7-0  n  *-i  ;-i 


<  — ~  <Z)22(yj)-2)siny1.a(^-r))  + 

2  i-i  4 


(95) 


11  r 

E  ^  r 

j  -2 


sinaa(tt  -r) 

c? 


da  ]  T«4)  +  I  y^(r)  v/l-r2  }. 
4  j 


The  treatment  of  the  integrals  in  Eqns.  (79)  and  (80)  with  kernels  k]2(r,s)  and 
k2I(r,s)  in  converting  them  to  summations  is  very  similar  to  those  with  kernels  ku(r,s) 
and  k2/r,s)  obtained  earlier.  The  difference  is  that  instead  of  having  to  deal  with  a 
step  function  in  the  integrand  as  shown  in  Eqn.  (67),  there  is  now  a  logarithmic  term 
(see  Eqn.  (71)).  Utilizing  the  following  identity  [24]. 


f, 


Tj(s)  log  |$-r  | 
v/1  -s2 


ds  -  TXr),  r  <  1, 
J 


(96) 


The  second  integral  in  Eqn.  (79)  may  be  expressed  as  follows 

P  k12{r,s)  f2{s)a  ds  -  na  Y  Bi j  (—  Y  CC  (—2—  Y 

>-o  n  *.j  j.i  2  j.i 

NX  V'  •  f-  cosa a(tk-r) 
^(y.lcosy.a^-r))  +  Y  cn.i  - — - da  - 

« -2  jA  O' 


(97) 


1(1.  *  log \Aa  |  *  f"  ffg?-1.  da)  ]  TUt)  .  1  121  I, 

4  J 0  a  A  j 

Again,  all  relevant  variables  are  defined  earlier  in  this  section  except  for  y0  = 
0.57721566490,  the  Euler’s  constant,  which  is  mentioned  in  section  3.2.1.  The  integral 
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(98) 


m««,-r)i  cosa-l  ££a 

Jo  a 

is  evaluated  using  subroutine  available  in  the  IMSL  software  library.  The  subroutine 
is  based  on  auxiliary  functions  described  in  [23].  As  for  the  other  integral  in  Eqn. 
(97),  expressions  derived  in  Appendix  C  gives 
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E  ci.  J 

n-  2  JA 


-  cos  a  a(tk-r) 


a" 


da 
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E  C12,2»  j 

n-1  *'A 


■-  cosaa(£A-r) 


a 


2  n 


da  +  £  ci2,2i,-i  f 

n  *2  */A 


•-  cosa a(th-r) 
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2n-l 


da. 


where 


E  C12,2n  J 

n-1 


•-  cosa(fA-r) 
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2n 


da  - 


52  C12,2n 

n  ■! 


[  cosA  (tk  -  r)  52 
j-i 


(-lV*1(£fc-r)2t'-1)(2n-2j)! 
(2n-l)!  A2"-2-'*1  + 


(100) 


sinAft*  -r) 


n-1 

E 


j- i 


(-iy«4 -r)^'1(2n -2/-1)! 
(2n-l)!  A2n'2J 


(-1)"*1 


(£*  -r)2"'1 
(2n-l)! 


si(A(£4  -r))]. 
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E  C12.2n-1  J 

n-2 


cosa  (tk  -r) 


A  a2«-l 


da  ■ 

"  - 1  -r)2°'1)(2n -2j  - 1)! 


(101) 


E  ci2.2»-i  tcosAU*-r)  £ 

n-2  j*x  (2n -2)!  -A  ^ 

A/.  .V*  (-1V(<4 -r)2j~l(2n-2j -Z)\ 

sinA(#4-r)  >  I _  + 

U  (2n-2)!  A2”-2^1 
-r)2""2 

(-1)"  li — 1 -  Ci(A(tk-r))). 

(2n  -2)!  * 


Using  the  same  procedure,  the  integral  with  kernel  k21(r,s)  in  Eqn.  (80)  may  be 
expressed  as 

n  JV 


P  k21(r,s)  /i(s)a  ds  -  Jta  £  I—  E  lE  (  ^ 

>.o  n  *.i  i.i  2 

E  Z>n(yt)cosyta(tt-r))  +  £  c2*w  jf”  cosaa(<*  r) 

j  _  i  a  •'A  rv1 


i-1 

7 


da  - 


(102) 


-L(Yo  +  log  |Aa  |  +  fAo(‘*'r)'  da)  ]  T/tk)  *  1  It 2  }. 

4  Jo  a  4  j 


But  as  mentioned  earlier,  k2/r,s)  =  -k21(r,s),  Eqn.  (102)  therefore  need  not  be 
computed,  but  can  be  obtained  by  simply  changing  the  sign  of  kj/rJJ  which  has 
already  been  calculated. 

We  have  shown  that  the  sine  integral  si(A( tk-r)),  which  contains  a  step  function, 
is  to  be  separated  from  the  kernel  and  integrated  in  closed  form.  A  similar  procedure 
is  applied  to  the  cosine  integral  Ci(A(tk-r)),  which  has  a  logarithm  term  (see  Eqns.  (84) 
and  (96)).  These  steps  are  taken  because  the  numerical  procedure’s  inability  to 
integrate  properly  either  the  step  function  or  the  log  function.  Recognizing  the 
presence  of  additional  sine  integrals  in  Eqns.  (93)  and  (100),  as  well  as  cosine 
integrals  in  Eqns.  (92)  and  (101),  one  might  suspect  whether  it  is  also  necessary  to 
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separate  these  terms  which  contain  the  discontinuous  functions  mentioned  above.  A 
closer  look  at  Eqns.  (92),  (93)  and  (100),  (101),  however,  reveals  that  the  sine  and 
cosine  integrals  have  accompanying  power  functions.  The  presence  of  these  power 
functions  cause  the  discontinuity  in  the  step  function  and  the  weak  singularity  in  the 
logarithm  function  to  be  both  continuous.  Thus  we  can  proceed  with  the  Gaussian 
integration  and  get  a  good  approximation  without  having  to  separate  these  sine  and 
cosine  integrals.  In  Appendix  D,  we  shall  illustrate  this  fact  by  way  of  examples  to 
show  the  dampening  effect  of  the  power  function  on  the  step  function  as  well  as  the 
logarithm  function  and  the  usage  of  Gauss’  formula  to  evaluate  integrals  involving 
these  functions. 

We  normalize  the  single-valuedness  condition  of  Eqn.  (55)  so  that  the  limits  of 
integration  become  from  -1  to  1 

JT /;(«)  ds  -  0,  i«l,  2,  (103) 

Substitute  the  representation  of  the  density  functions  in  terms  of  summation  of 
Chebyshev  polynomials  and  the  unknown  coefficients 

p  T(<s) 

Us)  -  E  — — .  i  -  1,  2,  (104) 

j‘°  VT-s* 


into  Eqn.  (103),  and  noting  that  T0  (s)  -  1,  we  obtain 


t  b„  f,  d.  ■  0,  i  -  1,  2. 

>-°  J 1  -s2 


(105) 


vi-s2 


Since  the  Chebyshev  polynomial  of  the  first  kind  are  orthogonal  with  respect  to  the 

weight  function  - .  ,  we  conclude  that 

\/T -s2 
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Bi o  -  0,  i  -  1,  2. 


(106) 


The  unknowns  in  the  linear  system  thus  become,  Blj  and  B2j  ,  j  =  1  ••••  p,  a  total  of 
"2p"  unknowns  To  solve  this  linear  system,  we  need  to  know  "2p"  distinct  values  in 
the  forcing  functions  as  the  right  hand  side  of  the  linear  system.  Let 

r  -  cos  ~  m  -  1  •••  p.  (107) 

m  2  p 

rw  m  =  1  ••••  p  are  the  points  where  the  forcing  functions  p1  (r)  and  p2  (r)  are  to  be 
evaluated,  thus  forming  the  right  hand  side  vector.  Based  on  the  foregoing  discussion, 
we  can  write  the  system  of  singular  integral  equations  as  a  system  of  linear  equations 
as  follows 


K  +  l 


„  PMJ  -£  TjUt >* 

j. i  n  k-i 
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where 


49 


(110) 
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11 
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sinag(fA-rm) 
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(113) 


Eqns.  (108)  and  (109)  form  a  ”2p"  by  "2p"  linear  system. 

It  should  be  painted  out  that  the  choice  of  collocation  points  rm,m  =  1  — •  p  in 
Eqns.  (108)  and  (109)  is  arbitrary.  The  only  restriction  is  that  rm  *  tk,  which  is 
imposed  out  of  numerical  convenience.  Because  rm  =  tk  would  make  any  closed  form 
evaluation  of  integrals  that  involves  oscillatory  sine  and  cosine  functions  invalid.  On 
the  other  hand,  the  choice  of  collocation  points  affects  convergence  of  the  final 
solution.  According  to  [25],  it  is  best  to  choose  them  so  that  they  are  zeroes  of 
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Chebyshev  polynomials.  Since  tk,  k  =  1  ••••  n  are  also  zeroes  of  Chebyshev  polynomials 
by  definition,  we  find  that  rm  =  tk  can  simply  be  avoided  by  letting  n  be  even  and  p 
odd,  or  vice  versa. 
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3.3  Stress  intensity  factors 


The  stress  intensity  factors  for  our  problem  sire  defined  as  follows 


k-^-a)  -  lim  v/2 J-a-x)  o^ix, 0),  (114) 

x— »  -a 


Aj(a)  ■  lim  ^2(x-a)  c^ix, 0),  (115) 

x— »  a 


k2(-a)  ■  lim  v^2(-a  -x)  c^Oc.O),  (116) 

x— *  -a 


A2(a)  «  lim  \l2(x-a)  t,0),  (117) 

x— »  a 


In  order  to  investigate  the  stress  singularity  near  the  crack  tip,  we  need  to 
obtain  stress  expressions  for  our  problem.  We  recognize  that  while  the  stresses  inside 
the  crack  give  us  the  system  of  singular  integral  equations,  these  equations  also 
provide  us  the  stress  expressions  outside  of  the  crack  at  y  =  0,  provided  the  density 
functions  are  known.  Thus  from  Eqns.  (46)  we  can  write 


—  *  +  (s,0)  -  f°  ^—dt  +  (aku(x,t)f1(t)  dt  *  kl2{x,t)f2(t)  dt, 

2px  t-x 


xOc  +  l). <J  (x,0)  .  f°  l^lLdt*  \-kJxmt)  dt+  [°k22<xmt)  dt, 

t  -x  •>-» 


ut) 


(118) 


x<  -a,  x>a. 

We  expect  the  stresses  to  be  singular  at  the  crack  tips.  Furthermore  from  Eqns.  (114) 
through  (117),  we  know  that  any  singular  terms  in  the  stresses  weaker  than  the 
square  root  singularity,  as  well  as  regular  terms,  will  vanish  because  of  the 
expressions  in  the  radical  and  the  limiting  process.  Also  recall  that  the  terms 
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involving  Fredholm  kernels  are  bounded  everywhere;  hence  the  limiting  process  will 
serve  to  eliminate  them  and  we  need  only  concern  ourselves  with  the  singular  terms. 
We  may  express  the  leading  terms  of  normal  stress  at  the  beginning  of  the  crack  tip 
as 


lim  c  (x,0)  -  lim  — ^ —  f°  l^Ldt. 

*->  *-»  -o  7t(K  +  l)  t~X 


(119) 


Recall  that 


W) 


Fjt) 


r,  i  "  1,  2. 


(120) 


(t+a)v2(a-t)va 

As  will  be  discussed  in  Appendix  F  regarding  the  behavior  of  Cauchy  integral,  we  can 
express  it  as  follows  (see  Eqn.  (F.  18)) 


I: 


fx(t)  4>!*(-a)eT  <f^(a)e  T  -i 

dt  -  iz - ( x+a )  2  -  n - ( x-a )  T  +  P(x), 


a  t-x 


sinii 

2 


sin— 

2 


(121) 


7c  i  -a)  (x+a)  7  -  ^(a)  (x-a)  T]  +  P(x). 


Also  from  Appendix  F,  we  can  show  that 

i.  '  *  *  /torn 

Fj(t)  «  4>j(£)  (a-t)T  +  i  4>J«)  «+a)T  +  <t>3<^)  (t+a)J  (a-tj1,  ’ 

In  Eqns.  (121)  and  (122),  P(x)  is  bounded  everywhere  except  possibly  at  the  end  points 
-a,  a,  where  it  has  singularities  no  higher  than  the  square-root  singularity  and  <J>/,  4>/ 
and  <j)3*  are  all  bounded  functions.  Thus,  we  obtain 

Fj(-a)  «  4>i(-a)  y/2a.  (123) 

We  can  easily  see  that  is  singular  at  x  =  -a  from  Eqns.  (119)  and  (121).  At  the 
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crack  tips,  the  stresses  will  be  dominated  by  the  singular  term.  We  may  write  Eqn. 
(119),  with  the  help  of  Eqn.  (123),  as  follows: 


lim  o^Oc.O) 


■  lim 


2  pj  4>j(-a) 
*  +  1 


x-t  -a  K  +  1 


Ft(-a) 
\/2a(-a  -x) 


(124) 


Substituting  Eqn.  (124)  into  Eqn.  (114),  we  obtain  the  Mode  I  stress  intensity  factor 
at  the  crack  tip  x  =  -a  as 


k^-a) 


F^-a). 


(125) 


Similar  to  Eqn.  (123),  from  Eqn.  (122)  we  obtain  F^a)  *  i  ^(a)  \/2a  .  The 

dominant  term  of  at  the  other  end  of  the  crack  tip  may  then  be  written  as 


lim  Ow(jc,0) 
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lim  - 

x-¥  a 


2jl! 
K  +  l 


i  4>j(a) 
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*  lim 

*— ►  a 


2pj  Fj(a) 

K  +  ^  i/2a(x  -  a) 


(126) 


Substituting  Eqn.  (126)  into  Eqn.  (115),  we  obtain 

*,(a)  -  - — — —  FJa).  (127) 

(k  +  1)  \fa 


The  mode  II  stress  intensity  factors  at  either  end  of  the  crack  tip  are  obtained  by  the 
same  process.  They  are 
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2u. 

k2(-a)  -  - i —  "  i-a). 

(k  +  1)  \[a 


(128) 


2u, 

k2(a)  -  -  - L_  E2(a). 

(k  +  1)  va 


(129) 


These  stress  intensity  factors  may  be  conveniently  evaluated  from 


F,(s)  -  jX-  T/s),  *  -  1,  2, 

j- 1 


(130) 


where  ,  i  =  1,  2,  j  =  1  •—  p,  come  directly  out  of  the  solution  of  the  linear  system 
from  Eqns.  (108)  and  (109). 


3.4  Crack  opening  displacements 

As  a  result  of  discussions  carried  out  in  section  3.2.3,  the  crack  opening 
displacements  may  readily  be  computed.  Combining  Eqns.  (85)  and  (86),  we  obtain 
the  y  component  of  the  COD  as 


,  Uj. i<*>  r  2 

u2(x,0+)  -  0,0c, 0-)  -  T  -BXj - l-(£)  . 

,-1  j  \  a 


(131) 


It  is  not  difficult  to  show  that  the  x  component  of  the  crack  opening  displacement 
would  be 


u2(x,0+)  -  u^x.O-)  -  £  -B2j  - - 

j- i  j 


Uj.  i<£> 

_ £_  i-(£)2. 


(132) 


Since  the  stresses  have  a  square-root  singularity,  the  crack  opening  displacements  are 
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bounded  at  the  crack  tips  and  they  vanish  as  the  square-root  of  the  distance  measured 
from  the  crack  tip,  as  shown  in  Eqns.  (131)  and  (132). 


3.5  Strain  energy  release  rates 

The  strain  energy  release  rate,  is  defined  as  one  half  the  rate  of  change  of 
strain  energy  per  distance  of  crack  propagation 

G  -  ±—  (133) 

2  da 

The  strain  energy  release  rate  may  be  seen  as  the  force  tending  to  open  the  crack 
surface.  Its  evaluation  requires  only  a  knowledge  of  the  stresses  and  displacements 
near  the  crack  tip.  From  reference  [26],  we  may  express  the  strain  energy  release 
rate  solely  in  terms  of  the  stress  intensity  factor  and  some  material  constants.  The 
opening  mode  and  sliding  mode  strain  energy  release  rate  are  found  to  be 


Gi 


ft(K  + 1) 

8  Pi 
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JC(K  +  1)  ,2 

-  - - -  «2> 
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(134) 


respectively.  The  total  strain  energy  release  rate  becomes 
G  -  7t(^1)  (k{  *  kl).  (135) 

8pj 

Let  G0  be  the  strain  energy  release  rate  of  a  crack  of  length  2a  in  an  infinite 
elastic  homogeneous  medium  with  the  same  material  properties  as  those  of  material 
1  under  uniform  normal  stress  =  a0  at  infinity  We  then  obtain 
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(136) 


Go 


ji(k+  1) 
8^ 


k1„  is  tiie  Mode  I  stress  intensity  factor  under  the  configuration  mentioned,  and  from  [27] 
we  know  that  -  cQ\[a.  Therefore 


(137) 


3.6  Convergence 

Near  the  end  of  section  3.2.3,  we  have  successfully  converted  the  system  of 
singular  integral  equations  to  a  linear  system  of  the  size  "2p"  by  "2p"  (Eqns.  (108)  and 
(109)),  where  " p ",  is  the  number  of  terms  in  the  summation  that  make  up  the  density 
functions  which  are  the  unknowns  in  the  SIE’s.  "p”  is  also  the  number  of  collocation 
points  which  we  evaluate  the  forcing  function  to  form  the  right-hand-side  vector  of  the 
linear  system.  A  look  into  the  previously  mentioned  equations  reveals  that  the 
number  of  terms  in  the  Gauss-Chebyshev  integration  formula  "n"  is  something  we 
should  also  determine  when  solving  the  linear  system.  Furthermore,  we  shall  decide 
the  "A"  value  which  is  the  integration  cut-off  point  for  the  semi-infinite  integral  in 
Eqn.  (65).  Therefore,  there  are  still  these  three  parameters  that  need  to  be  fixed  for 
the  linear  system.  We  say  the  solution  to  the  system  of  the  singular  integral 
equations  has  converged  when  changing  any  of  the  three  parameters  does  not  affect 
the  numerical  value  of  the  solution  of  the  linear  system. 

Of  the  three  parameters  that  we  mentioned  earlier,  the  Gauss  terms  "n"  and 
the  collocation  points  "p"  are  both  integers,  the  integration  cut-off  point  "A",  however, 
is  a  real  number.  While  the  relative  magnitude  of  the  two  integer  parameter  are 
easier  to  measure  in  attempting  to  achieve  convergence  of  the  solution,  since  one  could 
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start  with  an  arbitrary  small  integer,  we  need  to  establish  the  magnitude  of  the  third 
parameter.  Recall  in  Eqns.  (64)  and  (70)  that  "A"  is  the  value  which  makes  the 
integral 


.  sina(f-x) 

I  [£$<*) -d,‘(a)]  or  da,  ij  «  1,  2,  (138) 

A  cosa (t  - x ) 


negligible.  A  look  at  Eqns.  (B.59)  through  (B.62)  shows  that  the  last  term  of  d^Yo), 

v11 

i,  j  =  1,  2  is  in  the  form  of  - 1 _  0(0)  .  Therefore,  even  though  we  do  not  know 

4096  a11 

the  exact  expressions  of  D^'fa),  since  di;7o)  is  the  first  11  terms  of  DtJ'(oO  when  it  is 
expanded  asymptotically,  the  leading  term  of  Dtf(aJ  -  dtJ'(ctJ  must  be  a  term  like 

yl2 

- 1 -  0(0)  .  This  gives  us  a  basis  of  determining  the  value  "A".  Assume  that 

8192a12 

when 


8192  A12 


(139) 


where  e  is  a  very  small  number,  we  can  ignore  Eqn.  (138).  We  can  easily  deduce  that 


A  "  ...  '  •  (140) 

v/8192  e 

Granted  that  the  initial  assumption  of  £  might  not  render  an  "A"  such  that  Eqn.  (138) 
is  close  to  zero,  just  as  the  initial  guess  of  the  other  two  parameters  "p"  and  ”n”  might 
not  render  the  solution  anywhere  near  convergence,  however,  it  provides  something 
that  we  can  improve  on. 

From  the  stand  point  of  solving  the  singular  integral  equation,  we  say  that  the 
solution  has  converged  when  changing  the  "A"  value  produces  the  same  numerical 
value,  provided  the  other  two  parameters  have  been  determined.  With  the  perspective 
of  numerically  computing  the  integral  in  Eqn.  (65),  however,  the  stabilization  of 
solution  takes  on  a  different  meaning.  As  will  be  described  in  detail  in  Appendix  E, 
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when  "A"  reaches  to  a  point  where  it  makes  the  value  of  some  component  terms  in 
Du’( cl)  to  reach  the  limit  of  floating  point  overflow  for  a  real  number,  the  ter  u  P  ,'(cl) 
-  d12'( cl)  is  identically  zero  from  the  numerical  computation  stand  point.  Any  further 
increase  in  "A"  will  only  serve  the  same  purpose.  As  a  matter  of  fact,  it  is  even 
possible  to  have  "A"  so  large,  that  in  evaluating  the  kernel  kn(r,s)  the  second  integral 
in  Eqn.  (62)  can  be  completely  ignored.  But  to  do  so  is  to  defeat  the  whole  purpose 
of  going  through  the  exercise  of  asymptotic  expansions  detailed  in  Appendix  B.  But 
as  we  will  show  in  Appendix  E,  for  the  majority  of  the  cases  considered  in  this  study, 
it  does  not  take  a  large  "A"  to  make  the  terms  in  Dj'(a)  to  reach  the  machine 
constant.  From  the  above  discussion  we  conclude  that  as  £  gets  smaller,  Eqn.  (138) 
will  diminish  and  the  solution  of  the  SIE’s  must  converge. 

We  start  a  pilot  case  by  taking  for  the  value  of  the  material  constants, 
Poisson’s  ratio,  v  =  0.3,  the  non-homogeneity  parameter  of  material  2,  y  =  -  3.0, 
geometry  constants  such  as  the  thickness  of  material  1  ,h2-  10  a,  that  of  material  2, 
h2  -  a  (a  is  half  the  crack  length).  A  uniform  unit  normal  stress  over  the  entire  crack 
length  acts  as  the  forcing  function.  To  save  computer  time  a  liberal  value  of  £  =  1.0 
x  10  11  is  used.  To  our  pleasant  surprise,  we  find  that  the  solution,  even  though  far 
from  being  converged  for  the  "p"  and  "n"  parameters  we  have  selected,  does  not 
depend  on  ”n",  the  number  of  terms  in  Gauss-Chebyshev  integration  formula. 
Figure  4  shows  one  of  the  parameters  obtained  as  part  of  the  solution,  the  Mode  I 
stress  intensity  factor,  changes  as  the  number  of  collocation  points  "p"  increases,  but 
remains  virtually  constant  with  respect  to  the  changes  in  the  number  of  Gauss  terms 
”n". 

Having  fixed  one  parameter,  namely  the  gauss  terms,  in  the  quest  of  a 
converged  solution  of  our  linear  system,  the  next  logical  step  would  be  to  see  how  the 
solution  changes  with  respect  to  variations  in  the  number  of  collocation  points. 
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Figure  5  shows  that  as  the  number  of  collocation  points  increases,  the  solution 
converges  quickly.  A  number  of  no  more  than  24  collocation  points  is  sufficient  to 
make  the  solution  to  reach  the  asymptote  as  shown  in  Figure  5.  The  final 
convergence  of  the  problem  is  reached  when  the  solution  remains  accurate  to  a  certain 
digit  with  respect  to  the  variance  of  the  last  parameter  e. 

To  ensure  convergence  of  the  solution  for  all  the  cases  computed  in  this  study, 
every  number  is  computed  at  least  six  times;  each  time  varying  one  parameter.  The 
parameters  used  for  each  computation  is  listed  in  Table  I.  An  accuracy  of  four 
significant  figures  is  maintained  for  each  of  the  six  computations  in  every  case 
calculated  in  this  work. 

As  we  have  mentioned  earlier  and  as  will  be  explained  in  Appendix  E  in  detail, 
for  most  cases  it  does  not  take  a  very  large  "A"  for  Eqn.  (138)  to  become  negligible. 
From  the  vantage  point  of  numerical  computation  we  want  to  keep  "A"  as  small  as 
possible,  provided  we  can  still  obtain  a  converged  solution.  The  reason  for  this  is  that 
a  smaller  "A"  means  less  effort  in  integrating  Eqn.  (89).  From  Eqn.  (140)  we  can 
easily  see  that  "A"  is  indeed  rather  small.  For  y  =  1.0  and  e  =  1.0  x  10'11  we  need  only 
to  integrate  Eqn.  (89)  from  0  to  A  =  3.90.  When  y  =  0.1,  A  becomes  0.39  according  to 
Eqn.  (140).  But  there  is  a  limit  below  which  the  "A"  value  will  make  Eqn.  (138)  no 
longer  negligible.  Figure  6  illustrates  this  by  showing  a  plot  of  the  normalized  strain 
energy  release  rate  vs.  the  nonhomogeneity  constant  (y)  for  the  case  hj/a  =  100,  h2/a 
=  0.5,  v  =  0.3.  We  see  a  pronounced  dip  in  the  computed  parameter  when  |y|  <  1.0, 
where  it  should  have  been  a  very  smooth  transition  as  y  varies. 

Let  us  keep  in  mind  that  this  deviation  occurs  because  we  want  to  keep  "A" 
small  to  reduce  numerical  computing  effort.  If  we  have  kept  the  same  "A"  (the 
largest)  for  the  worst  case,  namely  |y|  =  3.0,  and  use  it  for  computing  any  other  case 
which  has  a  smaller  y,  the  deviation  would  not  have  happened.  We  would  at  most 
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wind  up  with  a  less  efficient  use  of  computing  time.  But  to  put  the  anomaly  just 
described  in  a  different  perspective,  let  us  explain  it  as  follows. 

As  |y|  gets  smaller,  the  magnitude  of  every  term  in  the  Fredholm  kernel  is 
getting  smaller  at  the  same  time.  But  the  basis  we  use  to  determine  "A",  is  built  on 
dropping  certain  terms  when  the  leading  term  is  less  than  an  absolute  number  e,  as 
shown  in  Eqn.  (139).  Therefore  the  relative  importance  of  these  terms  that  are  being 
dropped  increases  as  |y|  gets  smaller  and  smaller,  hence  the  deviations  as  shown  in 
Figure  6.  These  deviations  are  something  on  which  changing  the  parameters  as 
shown  in  Table  I  has  no  affect;  meaning  we  will  still  have  convergence  on  the  surface. 
Only  when  we  examine  any  of  the  parameters  that  we  seek,  such  as  SIF  or  strain 
energy  release  rate,  with  respect  to  the  changes  in  the  nonhomogeneity  constant,  do 
we  notice  the  aberrations. 

The  remedy  for  this  situation  is  of  course  trivial;  we  could  simply  use  a  larger 
"A".  For  those  cases  of  nonhomogeneity  constant  that  fall  into  this  range,  we  merely 
substitute  a  y  value  in  Eqn.  (140)  that  is  outside  of  the  range,  |y|  =  1.5,  to  be 
conservative  in  this  particular  example,  instead  of  using  the  actual  y  value. 
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Table  I.  Different  parameters  used  for  each  computation  in  order  to  obtain  a 
converged  solution. 


No. 

e 

P 

n 

1 

1.0  x  1011 

22 

51 

2 

1.0  x  1011 

22 

61 

3 

1.0  x  10  n 

24 

51 

4 

1.0  x  10  u 

24 

61 

5 

1.0  x  10-13  (1) 

24 

51 

6 

1.0  x  10'“  (1) 

24 

61 

Note:  (1)  Except  for  the  two  cases  (a).  h}  =  100a,  h2  =  0.5a  (b)  h2  =  100a,  h2  = 
0.25 a,  where  e  =  1.0  x  10'1S. 

(2)  These  parameters  are  used  for  all  the  cases  compiled  in  Chapter  4 
except  as  noted  in  (1). 

(3)  A  uniform  accuracy  of  four  significant  figures  convergence  is  obtained 
for  each  computation  with  the  parameters  used  as  lisi  i. 
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Figure  4.  Mode  I  SIF  vs.  number  of  Gauss  terms  with  collocation  points  2,  4,  6,  8 
for  e  =  1.0  x  10  u.  Material  constants:  v  =  0.3,  y  =  -3.0.  Geometry  constants:  h}  = 
10  a,  h2  =  a. 


2.46 


constants:  v  =  0.3,  y  =  -3.0.  Geometry  constants:  h}  =  10a,  h2  =  a. 
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number  of  collocation  points, 
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Figure  6.  Example  to  show  that  too  small  a  value  of  "A"  could  result  in  Eqn.  (138) 
not  negligible. 
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nonhomogeneity  constant 


3.7  Direction  of  crack  growth 

We  observe  that,  due  to  lack  of  symmetry  in  our  problem  since  material  2  is 
non-homogeneous  in  the  y-direction,  and  mixed-mode  condition  exists  regardless  of 
whether  the  loading  is  of  one  mode  or  mixed.  That  is  why  we’ll  see  in  Chapter  4  that 
both  kj  and  k2  are  non-zero  for  all  the  loading  conditions  computed.  With  this 
mixed-mode  conditions,  it  is  apparent  that  the  crack  will  not  advance  along  the 
interface.  In  order  to  predict  the  direction  in  which  the  crack  grows,  one  needs  to 
examine  the  stress  state  in  the  near  field  of  crack  tip.  The  stress  state,  in  polar 
coordinates  (see  Figure  7),  in  the  neighborhood  of  the  crack  tip  can  be  written  as 
follows  [28]: 

o_  -  — 1— cos— [A.(l  +sin2— )  +  —  A, sinG -2^, tan +0(r  2), 

VST  2  22  2 


CTee 


— ]L_cosi!.[A:1cos2^.--^.A2sin8]  +  0(rT), 

v/2T  2  2  2 


(141) 


«  — - — cos— [A1sin0-A2(3cos0-l)]  +  O(rT). 

2\f2r  2 

Following  [29],  [30],  the  direction  of  the  probable  crack  propagation  may  be 
determined  by  the  plane  of  the  maximum  cleavage  stress  ow(r,Q)  around  the  crack  tip. 
We  obtain  this  direction  either  by  solving  the  following  equation 

J_[v^7 o„(r,0)]  -  0.  (142) 

00 

which  is  to  find  the  direction  where  aw(r,Q)  remains  stationary,  or,  if  we  recognize 
that  it  is  equivalent  to  finding  the  directions  of  the  principal  stress,  simply  solving 
Oref r,  Q)  =  0.  Whichever  case  is  used,  we  obtain 
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(143) 


cos-ifAtSine+iJSsinO-l)]  *  0. 

2 

0 

cos_  -  0  gives  0  -  ±k,  which  corresponds  to  the  surface  of  the  crack.  The  second 
2 

term  in  Eqn.  (143)  prc  vides  the  non-trivial  solution  and  it  may  be  expressed  in  the 
following  quadratic  form  so  that  it  is  more  easily  solved 

tan2—  -_^L  tan—  -  —  •  0,  (144) 

2  2k2  2  2 

Note  that  the  direction  of  probable  crack  growth  is  dependent  only  on  the  stress  state 

i 

around  the  crack  tip.  Therefore  and  subsequent  terms  need  not  be  included  in 
the  foregoing  analysis  because  they  vanish  around  the  crack  tip. 


67 


Figure  8.  Probable  direction  of  crack  propagation. 


Chapter  4 

Results  And  Discussions 


4.1  Introduction 

In  applying  the  solution  of  our  problem  to  actual  cases,  we  consider  two 
material  combinations.  One  that  has  a  non-homogeneous  material  which  is  "softer" 
than  the  homogeneous  material,  as  can  be  represented  by  y  <  0.  For  example,  a  thin 
film  of  gold,  silver  or  lead  to  be  deposited  on  Nickel  or  steel  bearing  and  sliding  parts 
as  dry  film  lubricant[l][2]  has  this  kind  of  characteristics.  Gold  has  a  modulus  of 
elasticity  of  10.8  x  106  psi.,  Nickel  and  steel  have  moduli  of  elasticity  of  32  x  106  psi. 
and  29  x  106  psi.  respectively.  On  the  other  hand,  refractory  heat-shield  materials  on 
metal  substrates  may  or  may  not  have  a  surface  Young’s  moduli  higher  than  those  of 
the  substrate.  Refractory  materials  have  Young’s  moduli  anywhere  from  10  to  70  x 
106  psi.  The  range  of  Young’s  modulus  for  the  metal  substrate  could  be  even  more 
diverse.  Thus  we  may  have  a  homogeneous  material  either  stiffer  or  softer  than  the 
nonhomogeneous  layer,  depending  on  the  particular  application. 

The  material  and  geometry  constants  used  in  the  computation  of  the  cases  in 
this  chapter  are  designed  to  cover  a  broad  range  of  possible  material  and  geometry 
combinations.  The  non-homogeneity  constant  used  in  the  computation,  ranges  from 
-3  to  3.  As  for  Poisson’s  ratio,  a  value  of  0.3  was  used  in  the  computation.  A  separate 
investigation,  however,  on  the  effect  of  the  changes  in  Poisson’s  ratio  is  done  on  one 
of  the  geometry  combinations.  The  different  cases  of  geometry  and  material  constant 
combinations  computed  is  shown  in  Table  II.  Only  two  loadings  cases  were 
investigated,  namely  the  uniform  normal  stress  and  uniform  shear  applying  at  the 
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Table  II.  Material  and  geometry  constants  combinations  computed  in  the  present 
work. 


Case  No. 

hj/a 

hj/a 

Poisson’s  ratio  v 

1 

100 

2 

10 

3 

100 

2 

4 

1 

5 

0.5 

6 

0.25 

7 

10 

0.3 

8 

4 

1 

9 

2 

10 

1 

11 

100 

1 

0.01,  0.1,  0.2,  0.3,  0.4, 

0.499 

top  and  bottom  surfaces  of  the  crack  (Figure  9).  Due  to  the  myriad  of  loading 
combinations  conceivable,  these  two  loadings  computed  illustrate  perhaps  the  most 
general  cases.  In  addition,  the  result  of  these  two  loading  cases  for  when  the  two 
materials  are  very  thick  provided  verification  and  comparison  of  results  of  a  similar 
problem  with  two  half  planes. 

4.2  Results  and  discussions 

The  first  thing  to  do  after  obtaining  the  solution  of  the  system  of  integral 
equations  is  to  compute  a  case  with  parameters  very  close  to  one  whose  solution  has 
been  known  and  to  compare  the  results.  Reference  [16]  solved  an  interface  crack 
problem  of  a  similar  nature  to  what  is  done  in  the  present  work.  The  only  differences 
are  that  in  [16]  the  two  materials  are  both  half-spaces  rather  than  of  finite 
thicknesses  and  the  Poisson’s  ratio  exhibits  the  same  nonhomogeneous  properties  as 
that  of  the  Young’s  modulus  for  material  2,  whereas  in  this  study  the  Poisson’s  ratio 
is  treated  as  a  constant.  To  make  a  proper  comparison,  we  can  make  the  thickness 
of  both  materials  to  be  very  large,  e.g.,  h2  =  h2  -  100a.  As  for  Poisson’s  ratio,  a  value 
of  0.3  is  used.  We  will  show  later  in  this  chapter  that  the  effect  of  the  changes  in 
Poisson’s  ratio  is  not  significant. 

Figure  10  and  Figure  11  show  the  Mode  I  and  II  stress  intensity  factors  under 
the  loading  of  normalized  uniform  normal  stresses  of  the  above  mentioned  two  cases. 
Figure  12  and  Figure  13  show  the  same  stresses  intensity  factors  under  normalized 
shear  stresses.  Comparison  of  the  these  results  shows  that  they  can  be  very  closely 
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Figure  9.  Loadings  used  in  the  present  work,  (1).  uniform  normal  stress, 
(2).  uniform  shear. 
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Figure  10.  Comparison  of  Mode  I  SIF  for  interface  crack  between  (1).  two  infinite  half 
planes,  and  (2).  hi=h2=100a  under  loading  of  uniform  normal  stress. 
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Figure  11.  Comparison  of  Mode  II  SIF  for  interface  crack  between  (1).  two  infinite 
half  planes,  and  (2).  h1=h2=100a  under  loading  of  uniform  normal  stress. 
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Figure  12.  Comparison  of  Mode  I  SIF  for  interface  crack  between  (1).  two  infinite  half 
planes,  and  (2).  h^h^lOOa  under  loading  of  uniform  shear. 
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Figure  13.  Comparison  of  Mode  II  SIF  for  interface  crack  between  (1).  two  infinite 
half  planes,  and  (2).  h^h^lOOa  under  loading  of  uniform  shear. 
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correlated  with  differences  perhaps  attributed  to  the  difference  in  Poisson’s  ratio. 

We  next  examine  the  various  cases  computed  in  this  work.  Figure  14  through 
Figure  25  show  the  normalized  strain  energy  release  rate  and  Mode  I  and  Mode  II 
stress  intensity  factors  for  the  cases  where  h2  =  100a  remains  unchanged,  only  h2 
varies,  under  the  two  loadings  described  in  Figure  9. 

First  of  all,  all  the  parameters  computed  for  cases  (1),  h2  =  100a,  h2  =  100a,  and 
(2).  hj  =  100a,  h2  =  10a,  are  identical  to  the  third  digit  on  the  right  hand  side  of  the 
decimal  point.  This  is  why  results  of  case  (1)  are  not  shown  in  Figure  14  through 
Figure  19.  And  it  tells  something  about  the  nature  of  the  cases  when  both  materials 
are  thick;  as  the  nonhomogeneous  material  gets  thicker,  beyond  h2  >  10  (may  be  even 
less  than  10),  any  further  increase  in  h2  has  no  effect  on  the  outcome  of  the  solution. 

When  y  <  0,  it  is  as  if  any  material  beyond  h2  >  10  does  not  exit,  because  large  h2  and 
the  negative  sign  in  the  exponential  have  made  Pj  close  to  zero.  On  the  other  hand, 
the  positive  y,  together  with  large  h2  has  made  materials  beyond  h2>  10  to  be  like  a 
rigid  body.  It  explains  why  we  see  in  Figure  14  through  Figure  19  that  the  all 
parameters  tend  to  converge  at  the  two  extremes  of  the  nonhomogeneity  constant 
where  they  are  essentially  the  same  as  the  half  planes  case.  Only  when  y  is  small  do 
we  see  the  deviation  from  the  half  planes  case. 

The  Mode  I  stress  intensity  factor  under  both  loadings  exhibits  a  downward 
trend  as  the  nonhomogeneity  constant  y  increases  from  negative  to  positive.  This  can 
be  explained  partly  with  the  help  of  the  discussions  on  the  relative  stiffness  of  the 
nonhomogeneous  material  as  its  thickness  increases.  The  changes  in  the 
nonhomogeneity  constant  is  showing  a  similar  effect.  Large  y  on  the  negative  side 
means  the  crack  is  closer  to  the  free  surface;  smaller  ligament  thickness  indicates 
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easier  tear  along  the  interface,  therefore  a  larger  mode  I  stress  intensity  factor  as 
shown  in  Figure  15  and  Figure  18.  This  argument  explains  clearly  the  upward  trend 
of  kj  in  the  same  figures  as  the  nonhomogeneous  layer  gets  thinner.  Since  case  (2' 
hj  =  100a,  h2  =  10a,  is  identical  to  the  half  planes  case  as  far  as  the  number  of  digits 
in  the  parameters  we  obtain,  it  is  only  appropriate  that  the  stress  intensity  factors 
exhibit  the  same  traits  as  the  half  planes  case.  That  is  why  for  y  =  0  we  see  for  case 
(2).  the  Mode  I  SIF,  k/a)  =  1.0  and  Mode  II  SIF,  k/a)  =  0  under  uniform  normal 
stress,  as  shown  in  Figure  15  and  Figure  16.  Similarly  for  the  same  case  (2).,  Mode 

I  SIF,  k/a)  =  0  and  Mode  II  SIF,  k2(a)  =  1.0  when  the  loading  is  uniform  shear,  as 
shown  in  Figure  18  and  Figure  19. 

Due  to  the  nonhomogeneity  of  the  medium,  the  stress  intensity  factors  exhibit 
mixed  mode  condition  even  though  the  loadings  are  of  one  mode.  But  in  general, 
Mode  II  SIF  is  of  secondary  importance  as  compared  to  Mode  I  under  uniform  normal 
stress,  just  as  Mode  I  SIF  is  secondary  as  compared  to  Mode  II  when  the  loading  is 
uniform  shear.  These  facts  can  be  verified  by  the  relative  magnitude  of  k2  and  k2 
under  respective  loading  as  shown  in  Figure  15  through  Figure  19.  Furthermore, 
because  the  geometry  of  the  problem  is  symmetric  with  respect  to  the  y-axis,  the  shear 
stresses  are  anti-symmetric  under  uniform  normal  stress  loading.  Therefore  the  Mode 

II  SIF  under  uniform  normal  stress  has  different  signs  at  each  end  of  the  crack  tip. 
This  can  also  help  to  explain  why  k/a)  goes  through  a  sign  change  as  y  changes  from 
negative  to  positive  in  Figure  16.  We  know  that  when  y  <  0,  the  homogeneous 
material  is  stiffer  than  the  nonhomogeneous  material.  Whereas  the  reverse  is  true 
when  y  >  0.  We  can  envision  the  change  of  sign  in  shear  stress  as  y  goes  from 
negative  to  positive  as  if  the  geometry  of  the  problem  has  been  rotated  180  degrees. 
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Because  what  used  to  be  stiffer  (material  1)  when  y  <  0  is  "softer"  when  y  >  0,  the  sign 
of  the  shear  stresses  at  x  =  -a  when  y  <  0,  which  is  positive,  is  exactly  what  it  should 
be  at  x  =  a  when  y  >  0.  But  as  h2  gets  thinner  and  thinner,  the  tendency  of  this  Mode 
II  SIF  under  normal  stress  to  change  sign  as  y  goes  from  negative  to  positive  because 
of  the  relative  material  stiffness  explained  earlier  as  shown  in  Figure  16  is  being 
offset  by  the  less  and  less  ligament  between  the  crack  and  the  free  surface.  As  a 
result  the  tendency  to  change  sign  is  less  as  h2  becomes  thinner.  This  is  reflected 
clearly  in  Figure  22  where  we  see  for  h2  =  a,  k2  changes  its  sign  at  y  =  2,  whereas  for 
h2  =  0.5a  and  h2  =  0.25a,  k2  remains  negative. 

It  should  be  pointed  out  that  for  uniform  shear  loading,  is  symmetric  with 
respect  to  the  y-axis  and  ow  becomes  anti-symmetric.  As  a  result  k2  is  positive  at 
either  tip  of  the  crack  whereas  k1  undergoes  a  sign  change  from  one  crack  tip  to 
another.  But  examination  of  the  vertical  crack  opening  displacement  under  this 
loading  (see  Figure  41)  reveals  that  ( v *  -  V)  is  negative  on  the  left  half  of  the  crack. 
This  means  the  upper  crack  surface  has  penetrated  into  the  lower  crack  surface, 
which  is  not  permissible.  What  actually  would  have  happened  is  that  the  crack  would 
have  opened  for  only  a  portion  of  the  crack  length,  and  would  remain  closed  for  the 
remainder.  The  problem  needs  to  be  reformulated  if  we  want  to  solve  it  under  this 
loading. 

The  normalized  strain  energy  release  rate  would  reflect  the  combined  effect  of 
both  kj  and  k2  in  each  loading  for  it  is  a  function  of  (k2  +  k2  )  as  discussed  in  Chapter 
3. 
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Figure  14.  Normalized  strain  energy  release  rate  of  interface  crack  for  (1).  hi=100a, 
h2=10a,  (2).  h3=100a,  h2=2a,  under  loading  of  uniform  normal  stress. 
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Figure  15.  Normalized  Mode  I  SIF  of  interface  crack  for  (1).  hi=100a,  h2=10a,  (2). 
hj=100a,  hz=2a,  under  loading  of  uniform  normal  stress. 
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Figure  16.  Normalized  Mode  II  SIF  of  interface  crack  for  (1).  hj=100a,  h2=10a,  (2). 


hi=100a,  h2=2a,  under  loading  of  uniform  normal  stress. 
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h^lOOa,  hz=2a,  under  loading  of  uniform  shear. 
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Loading  =  uniform  shear  stress 


Figure  19.  Normalized  Mode  II  SIF  of  interface  crack  for  (1).  h^lOOa,  h2=10a,  (2). 
hjslOOa,  h2=2a,  under  loading  of  uniform  shear. 
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Figure  20.  Normalized  strain  energy  release  rate  of  interface  crack  for  (1).  h^lOOa, 
h2=a,  (2).  hi = 100a,  h2=0.5a,  (3).  h1=100a,  h2=0.25a,  under  loading  of  uniform  normal 
stress. 
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Figure  21.  Normalized  Mode  I  SIF  of  interface  crack  for  (1).  hj=100a,  h2=a,  (2). 


h^lOOa,  h2=0.5a,  (3).  h1=100a,  h2=0.25a,  under  loading  of  uniform  normal  stress. 
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Figure  22.  Normalized  Mode  II  SIF  of  interface  crack  for  (1).  hj=100a,  h2=a,  (2). 
h^lOOa,  h2=0.5a,  (3).  h^lOOa,  h2=0.25a,  under  loading  of  uniform  normal  stress. 
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Figure  23.  Normalized  strain  energy  release  rate  of  interface  crack  for  (1).  h^lOOa, 
h2=a,  (2).  hjslOOa,  h2=0.5a,  (3).  h1=100a,  1^=0.258,  under  loading  of  uniform  shear. 
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Figure  24.  Normalized  Mode  I  SIF  of  interface  crack  for  (1).  hj=100a,  h2=a,  (2). 
h^lOOa,  h2=0.5a,  (3).  h^lOOa,  h2=0.25a,  under  loading  of  uniform  shear. 
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Figure  25.  Normalized  Mode  II  SIF  of  interface  crack  for  (1).  h^lOOa,  h2=a,  (2). 
hjslOOa,  h2=0.5a,  (3).  h^lOOa,  h2=0.25a,  under  loading  of  uniform  shear. 
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Figure  26  through  Figure  31  show  the  normalized  strain  energy  release  rate 
and  Mode  I  and  Mode  II  stress  intensity  factors  when  h2  =  a  and  h2  varies  under  the 
two  loading  conditions.  Again,  the  Mode  I  stress  intensity  factor  under  uniform 
normal  stress  decreases  as  y  goes  from  negative  to  positive.  It  also  increases  with 
decreasing  ligament  size  in  material  1,  as  can  be  seen  in  Figure  27.  The  Mode  II  SIF 
under  uniform  normal  stress  shows  the  same  tendency  to  change  sign  as  y  goes  from 
negative  to  positive.  For  the  same  reason  in  the  cases  where  h2  =  100a  and  h2  varies, 
Figure  28  shows  that  the  tendency  for  the  sign  change  is  discouraged  as  h2  becomes 
thinner.  It  is  worth  noting  that  only  when  h1  -h2-  a  and  y  =  0  do  we  see  one  mode 
as  signified  by  k2  =  0  under  uniform  normal  stress  in  Figure  28  and  k2  =  0  under 
uniform  shear  as  seen  in  Figure  30.  The  rest  of  the  cases,  either  because 
nonsymmetric  geometry  or  nonsymmetric  material  constant  (nonhomogeneity),  are  of 
mixed  mode. 

The  effect  of  the  changes  in  Poisson’s  ratio  on  the  strain  energy  release  rate 
and  stress  intensity  factors  for  the  case  h2  =  100a,  h2  =  a  under  the  two  loadings  are 
shown  in  Figure  32  through  Figure  37.  We  note  that  for  plane  strain  k  =  3  -  4v, 
therefore  k  is  equal  to  1  and  3  respectively  for  Poisson’s  ratio  of  0.5  and  0.  Either  of 
these  two  values  of  k  causes  some  term  in  the  denominator  to  become  zero  during  the 
course  of  numerical  computation  (see  Eqn.  (1)  in  Chapter  2).  To  avoid  dividing  by 
zero,  Poisson’s  ratios  of  0.01  and  0.499  are  used  instead.  From  Figure  32  through 
Figure  37,  we  see  that  changes  in  Poisson’s  ratio  does  not  have  a  pronounced  effect 
on  the  parameters  calculated,  except  for  very  large  values  of  y,  as  evidenced  by  the 
closeness  of  the  curves.  These  comparisons  also  show  that  when  the  materials  are 
homogeneous,  the  stress  intensity  factors  are  independent  of  the  Poisson’s  ratio.  It 
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Figure  26.  Normalized  strain  energy  release  rate  of  interface  crack  for  (1).  h^lOa, 
h2=a,  (2).  h!=4a,  h2=a,  (3).  h1=2a>  h2=a,  (4).  hx=h2=a,  under  loading  of  uniform  normal 
stress. 
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h2=a,  (3).  hj=2a,  h2=a,  (4).  hi=h2=a,  under  loading  of  uniform  normal  stress. 
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Figure  29.  Normalized  strain  energy  release  rate  of  interface  crack  for  (1).  hj=10a, 
h2=a,  (2).  hx=4a,  h2=a,  (3).  hx=2a,  h2=a,  (4).  hx=h2=a,  under  loading  of  uniform  shear. 
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Figure  30.  Normalized  Mode  I  SIF  of  interface  crack  for  (1).  hjslOa,  h2=a,  (2).  h,=4a, 
h2=a,  (3).  h,=2a,  h2=a,  (4).  h^hj^a,  under  loading  of  uniform  shear. 
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Figure  31.  Normalized  Mode  II  SIF  of  interface  crack  for  (1).  hi=10a,  h2=a,  (2).  hj=4a, 


h2=a,  (3).  hj=2a,  h2=a,  (4).  h1=h2=a,  under  loading  of  uniform  shear. 
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is  evidenced  by  the  SIF’s  for  all  values  of  Poisson’s  ratio  converging  at  y  =  0.  Before 
we  attempt  to  justify  the  behavior  of  the  parameters  shown  in  Figure  32  th*  ^ugh 
Figure  37,  we  recognize  that  v  =  0  signifies  no  lateral  strain  and  v  =  0.5  represents 
no  volume  change,  which  is  also  the  maximum  lateral  strain  allowed.  Poisson’s  ratio 
of  value  larger  than  0.5  is  unlikely  because  in  the  case  of  tensile  loading  it  means 
there  is  a  volume  decrease.  Based  on  the  above  discussions,  when  y  <  0,  we  can 
interpret  v  =  0  to  be  stiffer  than  v  =  0.5,  therefore  a  smaller  SIF  as  seen  in  Figure  33 
and  Figure  37.  But  the  same  can  not  be  said  for  cases  when  y  >  0.  Actually  we  see 
just  the  reverse.  The  above  argument  can  also  be  applied  to  explain  the  small 
difference  between  the  parameters  computed  in  [16]  and  case  (1).  in  this  work.  As 
mentioned  earlier,  the  Poisson’s  ratio  is  treated  as  a  constant  in  this  work  and  varies 
the  same  way  as  does  the  shear  modulus  for  material  2  in  reference  [16].  As  a  result, 
everything  else  being  equal,  for  y  <  0  the  nonhomogeneous  material  in  [16]  is  actually 
stiffer  than  when  v  remains  constant  as  y  increases.  Therefore  Figure  10  and 
Figure  13  show  the  result  in  this  work  to  have  a  larger  stress  intensity  factor. 

Figure  38  through  Figure  41  show  crack  opening  displacements  for  the  case  h} 
=  100a,  h2  =  a  and  y  =  -3.0  whereas  Figure  42  through  Figure  45  show  COD’s  for  the 
same  geometry  except  y  =  3.0.  Figure  38  and  Figure  39  are  typical  of  the  crack 
opening  displacements  one  would  expect  under  uniform  normal  stress.  The  horizontal 
component  displays  a  movement  toward  the  center  of  the  crack  at  which  point  the 
horizontal  displacement  is  zero.  The  vertical  component  of  the  COD,  on  the  other 
hand,  has  an  elliptical  shape.  The  crack  opening  displacements  under  uniform  shear, 
however,  is  a  reverse  of  what  have  just  been  described.  The  horizontal  component  has 
an  elliptical  distribution  whereas  the  vertical  component  shows  a  change  in  direction 
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Figure  32.  Normalized  strain  energy  release  rate  of  interface  crack  for  hi=100a,  h2 
at  v=0.01,  0.1,  0.2,  0.3,  0.4,  0.499  under  loading  of  uniform  normal  stress. 
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hl=100a,  h2=a 

loading  =  uniform  normal  stress 


Figure  33.  Normalized  Mode  I SIF  of  interface  crack  for  hjslOOa,  h^a  at  v=0.01, 0.1, 
0.2,  0.3,  0.4,  0.499  under  loading  of  uniform  normal  stress. 
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Figure  34.  Normalized  Mode  II  SIF  of  interface  crack  for  h^lOOa,  h2=a  at  v=0.01, 0.1, 
0.2,  0.3,  0.4,  0.499  under  loading  of  uniform  normal  stress. 
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Figure  35.  Normalized  strain  energy  release  rate  of  interface  crack  for  h^lOOa,  hj=a 
at  v=0.01,  0.1,  0.2,  0.3,  0.4,  0.499  under  loading  oi  uniform  shear. 
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Figure  36.  Normalized  Mode  I  SIF  of  interface  crack  for  h^lOOa,  h.2=a  at  v=0.01, 0.1, 


0.2,  0.3,  0.4,  0.499  under  loading  of  uniform  shear. 
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0.2,  0.3,  0.4,  0.499  under  loading  of  uniform  shear. 
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of  the  relative  crack  surface  movements  on  either  half  of  the  crack,  as  seen  in 
Figure  40  and  Figure  41.  The  implication  of  this  behavior  is  quite  different  from  the 
horizontal  COD  under  uniform  normal  stress.  Negative  vertical  COD  is  not 
permissible  for  an  originally  closed  crack  because  it  means  the  crack  surfaces  have 
penetrated  each  other.  And  as  have  been  indicated  earlier,  in  this  case  the  crack 
would  remain  closed  for  some  portion  and  the  problem  would  have  to  be  reformulated. 

Figure  42  through  Figure  45  show  crack  opening  displacements  for  the  case 
otherwise  the  same  as  those  shown  in  Figure  38  through  Figure  41  except  y  =  3.0. 
One  easily  sees  the  difference  in  scale  between  corresponding  displacement  component 
in  these  two  groups.  The  former  represents  when  material  2  is  very  soft  (y  =  -3.0) 
therefore  crack  opening  displacements  are  much  larger  than  when  it  is  very  stiff  (y  = 
3.0).  There  is  a  peculiar  behavior  in  the  horizontal  component  of  COD  under  normal 
stress  as  well  as  the  vertical  component  od  COD  under  uniform  shear;  that  is  in 
addition  to  the  zero  movement  at  the  center  of  the  crack,  there  are  two  additional 
zeroes. 

The  angle  of  probable  crack  extension  at  crack  tip  x  =  aas  measured  from  the 
x  axis  as  a  function  of  the  nonhomogeneity  constant  is  shown  in  Figure  46  and 
Figure  47.  Because  the  uniform  shear  loading  case  shows  negative  Mode  I  SIF  which 
requires  problem  reformulation  due  to  portion  of  the  crack  not  being  open  as  discussed 
earlier,  only  the  uniform  normal  stress  loading  is  investigated.  Generally,  the  two 
figures  display  the  same  trend;  (1).  For  the  same  geometry,  a  shear  modulus  of  the 
nonhomogeneous  material  going  from  weak  to  stiff  will  see  the  angle  of  crack 
extension  change  from  large  to  small,  and  in  the  case  where  the  two  materials  of  the 


106 


'  i,  1  X  I  I - 1 - rr 

CO  CVj  T-  o  T-  CM  CO  Tt  ' 

OOOO  OOOO 

<111 


-n  -  +n 


Figure  38.  Crack  opening  displacements  (u*  -  u  )  of  interface  crack  for  hjslOOa,  hjsa, 
y=  -3.0  under  loading  of  uniform  normal  stress. 
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Figure  39.  Crack  opening  displacements  (v+  -  v  )  of  interface  crack  for  h^lOOa,  h^a, 
y=  -3.0  under  loading  of  uniform  normal  stress. 


108 


o 


Figure  40.  Crack  opening  displacements  (u*  -  u  )  of  interface  crack  for  hx=100a,  h^a, 
y=  -3.0  under  loading  of  uniform  shear. 
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Figure  41.  Crack  opening  1  ^placements  (v+  -  v')  of  interface  crack  for  h^lOOa,  hj=a. 


y=  -3.0  under  loading  of  uniform  shear. 
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Figure  42.  Crack  opening  displacements  (u*  -  u')  of  interface  crack  for  h^lOOa,  h^a, 
y=3.0  under  loading  of  uniform  normal  stress. 
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Figure  43.  Crack  opening  displacements  (v+  -  v')  of  interface  crack  for  h^lOOa,  hj=a, 
y=3.0  under  loading  of  uniform  normal  stress. 
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Figure  45.  Crack  opening  displacements  (v+  -  v')  of  interface  crack  for  hj^lOOa,  h^a, 
y=3.0  under  loading  of  uniform  shear. 
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same  thicknesses  from  positive  to  negative.  (2).  A  decrease  in  the  relative  ligament 
size  will  tend  to  discourage  the  trend  describe  in  (1).  Both  explain  the  fact  that  crack 
tend  to  grow  in  the  direction  where  the  combination  of  material  properties  and 
geometry  are  favorable  for  such  growth,  namely  a  less  stiff  material  property  or  a 
thinner  ligament  size.  Again,  the  case  h2  =  100a,  h2  =  100a  and  h2  =  100a,  h2  =  10a 
have  the  same  angle  of  extension  throughout  the  whole  range  of  nonhomogeneity 
constant  as  can  be  seen  in  Figure  46.  Similarly,  the  case  ht  =  100a,  h2  -  a  and  hx  - 
10a,  h2  =  a  show  exactly  the  same  behavior  as  seen  in  Figure  47. 

Table  III  through  Table  XX  list  the  Mode  I  and  Mode  II  stress  intensity  factors 
under  the  two  loadings  for  relative  dimensions  No.  1  through  No.  10  given  in  Table  II. 
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Figure  46.  The  angle  of  crack  extension  at  x=a  for  h^lOOa,  h2=100a,  10a,  2a,  a,  0.5a, 
0.25a,  respectively. 
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Table  HI.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hjslOOa,  bplOa,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kx(a)/ 

k^a)/ 

G/G0 

k^a)/ 

k^a)/ 

_  „ 

Poa 

P0a1/2 

_  1/2 
q0a 

„  _I/2 

9oa 

-0.01 

1.008 

1.011 

-0.002 

1.008 

0.002 

1.004 

-0.1 

1.087 

1.043 

-0.013 

1.033 

0.013 

1.016 

-0.25 

1.209 

1.099 

-0.034 

1.074 

0.032 

1.036 

-0.5 

1.437 

1.197 

-0.073 

1.138 

0.065 

1.065 

-0.75 

1.694 

1.296 

-0.116 

1.198 

0.097 

1.090 

-1.0 

1.978 

1.397 

-0.161 

1.255 

0.129 

1.113 

-1.25 

2.291 

1.499 

-0.210 

1.309 

0.159 

1.133 

m 

2.635 

1.602 

-0.261 

1.360 

0.188 

1.151 

-2.0 

3.422 

1.812 

-0.371 

1.460 

0.243 

1.184 

-2.5 

4.350 

2.027 

-0.489 

1.555 

0.293 

1.212 

-3.0 

5.432 

2.248 

-0.615 

1.648 

0.339 

1.238 
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Table  IV.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  bplOa,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0(1) 

k^a)/ 

k^a )/ 

G/G0(2) 

k^a)/ 

kjla)/ 

Po^ 

_  1/2 

Poa 

qoa1* 

„  „L/2 

9oa 

0.1 

0.950 

0.975 

0.012 

0.978 

-0.012 

0.989 

0.25 

0.870 

0.932 

0.028 

0.938 

-0.029 

0.968 

0.5 

0.775 

0.879 

0.050 

0.880 

-0.053 

0.937 

0.75 

0.715 

0.843 

0.067 

0.835 

-0.073 

0.911 

1.0 

0.674 

0.817 

0.081 

0.801 

-0.089 

0.891 

B 

0.624 

0.784 

0.103 

0.751 

-0.113 

0.859 

2.0 

0.596 

0.763 

0.120 

0.715 

-0.131 

0.836 

2.5 

0.577 

0.748 

0.132 

0.688 

-0.144 

0.817 

3.0 

0.564 

0.737 

0.142 

0.667 

-0.155 

0.802 
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Table  V.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hj=100a,  h2=2a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 


G/G0  kj(a)/  k2(a)/ 

„  „l/2  _  1/2 


1.361  1.166  -0.038 


1.415  1.189  -0.049 


1.513  1.228  -0.068 


1.699  1.299  -0.102 


1.913  1.376  -0.140 


2.159  1.458  -0.182 


2.439  1.545  -0.226 


2.754  1.637  -0.274 


3.495  1.831  -0.379 


4.393  2.037  -0.493 


5.457  2.253  -0.617 


uniform  shear 

G/G0 

kx(a)/ 

Ma)/ 

qoa" 

9oa 

1.091 

0.034 

1.044 

1.107 

0.043 

1.051 

1.134 

0.059 

1.063 

1.180 

0.085 

1.083 

1.227 

0.112 

1.102 

1.275 

0.140 

1.120 

1.322 

0.167 

1.138 

1.370 

0.193 

1.154 

1.464 

0.245 

1.185 

1.557 

0.294 

1.213 

1.649 

0.340 

1.238 

Table  VI.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  h2=2a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 


uniform  shear 


kj(a)/  ka(aV  G/G0  kj(a)/  k^a)/ 


1.298  1.139  -0.025  :  1.072  0.023  1.035 


1.220  1.105  -C.009  1.046  0.008  1.023 


.052  0.016 


0.037 


.965  0.056 


0.815  0.898  0.087 


0.732  0.849  0.110 


0.674  0.811  0.128 


0.633  0.783  0.141 


m 


-0.015  1.002 


0.965  -0.036  0.981 


0.927  -0.056  0.961 


0.860  -0.087  0.923 


0.803  -0.115  0.888 


0.754  -0.135  0.958 


0.714  -0.150  0.832 
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Table  VIE.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  lxjsa,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kj(a)/ 

k^a)/ 

G/G0 

k/a)/ 

Ma)/ 

Poa 

_ 

Poa 

_  „l/2 

9oa 

„  1/2 

9oa 

-0.01 

2.324 

1.513 

-0.186 

1.201 

0.134 

1.088 

-0.1 

2.389 

1.533 

-0.197 

1.214 

0.140 

1.093 

-0.25 

2.502 

1.567 

-0.216 

1.235 

0.152 

1.101 

-0.5 

2.706 

1.626 

-0.249 

1.271 

0.170 

1.114 

-0.75 

2.932 

1.689 

-0.284 

1.307 

0.189 

1.128 

-1.0 

3.180 

1.754 

-0.321 

1.345 

0.209 

1.141 

-1.25 

3.454 

1.823 

-0.360 

1.383 

0.228 

1.154 

mm 

3.755 

1.896 

-0.401 

1.422 

0.247 

1.167 

-2.0 

4.448 

2.051 

-0.491 

1.502 

0.286 

1.192 

-2.5 

5.273 

2.219 

-0.590 

1.584 

0.324 

1.216 

-3.0 

6.250 

2.400 

-0.698 

1.668 

0.361 

1.240 
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Table  VIII.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  h2=a,  v  -  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kt(a)/ 

k^a)/ 

G/G0 

kjCa)/ 

kj  (a)/ 

_ 

Poa 

_  _  1/2 

Poa 

q0a^ 

q0a^ 

0.1 

2.249 

1.490 

-0.173 

1.186 

0.126 

1.082 

0.25 

2.151 

1.458 

-0.156 

1.166 

0.115 

1.074 

0.5 

2.001 

1.409 

-0.129 

1.134 

0.097 

1.060 

0.75 

1.865 

1.362 

-0.103 

1.103 

0.079 

1.047 

1.0 

1.743 

1.318 

-0.079 

1.073 

0.062 

1.034 

1.534 

1.238 

-0.036 

1.016 

0.029 

1.008 

2.0 

1.365 

1.168 

0.001 

0.965 

-0.001 

0.982 

2.5 

1.226 

1.107 

0.034 

0.918 

-0.029 

0.958 

3.0 

1.113 

1.053 

0.063 

0.876 

-0.055 

0.934 
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Table  IX.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  h2=0.5a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

k/a)/ 

l^Ca)/ 

G/G0 

kj(a)/ 

k^a)/ 

_ 

Poa 

Poa1"2 

_  „l/2 

9oa 

qoa1* 

-0.01 

6.223 

2.400 

-0.682 

1.359 

0.319 

1.121 

-0.1 

6.333 

2.419 

-0.694 

1.369 

0.323 

1.125 

-0.25 

6.522 

2.452 

-0.715 

1.387 

0.330 

1.131 

-0.5 

6.852 

2.507 

-0.751 

1.417 

0.342 

1.140 

-0.75 

7.201 

2.565 

-0.788 

1.448 

0.354 

1.150 

l 

b 

7.572 

2.625 

-0.826 

1.480 

0.365 

1.160 

-1.25 

7.965 

2.686 

-0.866 

1.512 

0.377 

1.170 

m 

8.381 

2.749 

-0.906 

1.545 

0.389 

1.180 

o 

e4 

t 

9.288 

2.882 

-0.991 

1.612 

0.413 

1.201 

-2.5 

10.304 

3.022 

-1.082 

1.682 

0.437 

1.221 

-3.0 

11.441 

3.171 

-1.178 

1.753 

0.462 

1.241 

Table  X.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  h2=0.5a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

k^a)/ 

kjCa)/ 

G/G0 

kj(a)/ 

kj(a )/ 

Poa 

Poa”2 

q0a^ 

„  _V2 

9oa 

0.1 

6.091 

2.376 

-0.667 

1.346 

0.314 

1.117 

0.25 

5.917 

2.345 

-0.647 

1.329 

0.307 

1.111 

0.5 

5.640 

2.294 

-0.614 

1.301 

0.295 

1.102 

0.75 

5.380 

2.245 

-0.583 

1.273 

0.283 

1.092 

1.0 

5.133 

2.197 

-0.552 

1.247 

0.272 

1.083 

n 

4.682 

2.107 

-0.494 

1.196 

0.249 

1.065 

2.0 

4.281 

2.022 

-0.439 

1.148 

0.227 

1.047 

2.5 

3.924 

1.943 

-0.387 

1.103 

0.205 

1.030 

3.0 

3.605 

1.868 

-0.339 

1.062 

0.184 

1.014 
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Table  XI.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hi=100a,  h2=0.25a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

y 

G/G0 

k,(ay 

kjCa)/ 

G/G0 

kj(a)/ 

Ma)/ 

_ 

Poa 

_ 

Poa 

_  „l/2 

9oa 

q0al/2 

-0.01 

26.103 

4.642 

-2.134 

1.673 

0.540 

1.175 

-0.1 

26.369 

4.664 

-2.149 

1.683 

0.543 

1.178 

-0.25 

26.818 

4.699 

-2.176 

1.699 

0.548 

1.183 

-0.5 

27.586 

4.760 

-2.220 

1.725 

0.555 

1.190 

-0.75 

28.379 

4.822 

-2.264 

1.752 

0.563 

1.198 

-1.0 

29.297 

4.885 

-2.310 

1.779 

0.570 

1.206 

-1.25 

30.041 

4.949 

-2.356 

1.807 

0.578 

1.214 

g 

30.913 

5.014 

-2.403 

1.835 

0.586 

1.222 

-2.0 

32.740 

5.148 

-2.498 

1.893 

0.601 

1.238 

-2.5 

34.687 

.. 

5.286 

-2.597 

1.951 

0.616 

1.254 

-3.0 

36.760 

5.429 

-2.700 

2.012 

0.632 

1.270 

Table  XH.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOOa,  h2=0.25a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

y 

G/G0 

kjia)/ 

k^a)/ 

G/G0 

k/a)/ 

kj(a)/ 

Po*1* 

_ 

Poa 

„  „l/2 

9oa 

_  „l/2 

qoa 

0.1 

25.783 

4.616 

-2.115 

1.662 

0.537 

1.172 

0.25 

25.353 

4.581 

-2.090 

1.647 

0.532 

1.167 

0.5 

24.654 

4.523 

-2.048 

1.621 

0.525 

1.160 

0.75 

23.977 

4.467 

-2.006 

1.596 

0.518 

1.152 

1.0 

23.321 

4.411 

-1.966 

1.571 

0.510 

1.145 

n 

22.070 

4.303 

-1.886 

1.523 

0.495 

1.130 

2.0 

20.894 

4.198 

-1.809 

1.477 

0.481 

1.116 

2.5 

19.790 

4.097 

-1.734 

1.432 

0.466 

1.102 

3.0 

18.753 

3.999 

-1.662 

1.389 

0.452 

1.088 

Table  XIII.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^lOa,  h2=a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

k^a)/ 

k^a )/ 

G/G0 

k/a)/ 

kj(a)/ 

„  _  172 

Poa 

Poa 

9oa 

9oa 

-0.01 

2.325 

1.513 

-0.186 

1.205 

0.134 

1.090 

-0.1 

2.389 

1.533 

-0.197 

1.217 

0.141 

1.094 

-0.25 

2.502 

1.567 

-0.216 

1.238 

0.152 

1.102 

-0.5 

2.706 

1.626 

-0.249 

1.273 

0.171 

1.115 

-0.75 

2.932 

1.689 

-0.284 

1.310 

0.190 

1.129 

-1.0 

3.180 

1.754 

-0.321 

1.347 

0.209 

1.142 

-1.25 

3.454 

1.823 

-0.360 

1.385 

0.228 

1.155 

m 

3.756 

1.896 

-0.401 

1.424 

0.247 

1.168 

-2.0 

4.448 

2.051 

-0.491 

1.504 

0.286 

1.193 

-2.5 

5.273 

2.219 

-0.590 

1.586 

0.324 

1.217 

-3.0 

6.250 

2.401 

-0.698 

1.669 

0.361 

1.240 
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Table  XIV.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hi=10a,  hjsa,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

r 

G/G0 

kjCa)/ 

kj(a)/ 

G/G0 

k^a)/ 

k^a)/ 

_ 

Poa 

Poa 

_  \n 

9oa 

<lo*m 

0.1 

2.249 

1.490 

-0.173 

1.190 

0.126 

1.084 

0.25 

2.151 

1.458 

-0.156 

1.170 

0.115 

1.076 

0.5 

2.001 

1.409 

-0.129 

1.138 

0.097 

1.063 

0.75 

1.865 

1.362 

-0.103 

1.108 

0.079 

1.049 

1.0 

1.742 

1.318 

-0.079 

1.078 

0.062 

1.036 

Q 

1.535 

1.238 

-0.036 

1.023 

0.030 

1.011 

2.0 

1.365 

1.168 

-0.001 

0.973 

-0.001 

0.987 

2.5 

1.227 

1.107 

0.034 

0.928 

-0.030 

0.963 

3.0 

1.114 

1.053 

0.063 

0.886 

-0.056 

0.940 

Table  XV.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hj=4a,  kjsa,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

r 

G/G0 

kj(a)/ 

k^a)/ 

G/G0 

kj(a)/ 

k^a)/ 

P0aW 

„  1/2 

9oa 

„  „l/2 

9oa 

-0.01 

2.337 

1.518 

-0.184 

1.256 

0.135 

1.113 

-0.1 

2.401 

1.537 

-0.195 

1.267 

0.141 

1.117 

-0.25 

2.514 

1.571 

-0.214 

1.286 

0.153 

1.124 

-0.5 

2.717 

1.630 

-0.247 

1.319 

0.172 

1.136 

-0.75 

2.942 

1.692 

-0.282 

1.353 

0.191 

1.147 

-1.0 

3.191 

1.758 

-0.319 

1.388 

0.210 

1.159 

-1.25 

3.464 

1.826 

-0.358 

1.423 

0.230 

1.171 

Q 

3.765 

1.899 

-0.400 

1.460 

0.249 

1.182 

-2.0 

4.456 

2.053 

-0.490 

1.536 

0.288 

1.205 

-2.5 

5.281 

2.221 

-0.589 

1.614 

0.326 

1.228 

-3.0 

6.257 

2.402 

-0.697 

1.694 

0.363 

1.250 

Table  XVI.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h1=4a,  h^a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kj(a;/ 

kjCa)/ 

G/G0 

k^a)/ 

Ma)/ 

Po*1* 

_ 

Poa 

_  „l/2 

9oa 

q0al/2 

0.1 

2.262 

1.494 

-0.171 

1.242 

0.126 

1.107 

0.25 

-0.153 

1.224 

0.115 

1.100 

0.5 

2.014 

1.414 

-0.126 

1.195 

0.097 

1.089 

0  75 

1.880 

1.367 

-0.100 

1.167 

0.079 

1.077 

1.0 

1.758 

1.342 

-0.076 

1.140 

0.061 

1.066 

Q 

1.550 

1.245 

-0.032 

1.089 

0.027 

1.043 

2.0 

1.382 

1.175 

0.006 

1.042 

-0.005 

1.021 

2.5 

1.244 

1.115 

0.039 

0.999 

-0.035 

0.999 

3.0 

1.132 

1.062 

0.068 

0.957 

-0.062 

0.977 
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Table  XVH.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h!=2a,  h^a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

k^a)/ 

k^a)/ 

G/G0 

kx(a)/ 

kj(a)/ 

Po3^ 

Po3** 

q0a1/2 

Qoa 

-0.01 

2.478 

1.566 

-0.156 

1.455 

0.120 

1.200 

l 

O 

2.542 

1.585 

-0.168 

1.466 

0.127 

1.204 

-0.25 

2.653 

1.618 

-0.187 

1.483 

0.140 

1.210 

-0.5 

2.854 

1.675 

-0.221 

1.513 

0.161 

1.220 

-0.75 

3.077 

1.735 

-0.258 

1.544 

0.183 

1.229 

-1.0 

3.323 

1.799 

-0.296 

1.575 

0.204 

1.238 

-1.25 

3.594 

1.866 

-0.336 

1.607 

0.225 

1.247 

mm 

3.892 

1.936 

-0.379 

1.639 

0.246 

1.256 

-2.0 

4.577 

2.087 

-0.471 

1.706 

0.287 

1.274 

-2.5 

5.397 

2.252 

-0.571 

1.775 

0.328 

1.291 

-3.0 

6.367 

2.430 

-0.681 

1.846 

0.367 

1.308 
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Table  XV 111.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h!=2a,  h^a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kx(a)/ 

k2(a)/ 

G/G0 

k^a)/ 

Ma)/ 

_  1 
Poa 

Poa 

_  o1^ 

9oa 

q0aw 

0.1 

2.403 

1.544 

-0.143 

1.442 

0.111 

1.196 

0.25 

2.307 

1.514 

-0.125 

1.425 

0.098 

1.190 

0.5 

2.159 

1.466 

-0.096 

1.396 

0.077 

1.179 

0.75 

2.026 

1.422 

-0.069 

1.368 

0.057 

1.168 

1.0 

1.905 

1.380 

-0.044 

1.340 

0.037 

1.157 

n 

1.700 

1.304 

0.002 

1.284 

-0.002 

1.133 

2.0 

1.532 

1.237 

0.042 

1.229 

-0.038 

1.108 

2.5 

1.396 

1.179 

0.077 

1.174 

-0.071 

1.081 

3.0 

1.284 

1.128 

0.108 

1.120 

-0.101 

1.053 
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Table  XIX.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
h^a,  h^a,  v  =  0.3,  when  y  <  0. 


uniform  normal  stress 

uniform  shear 

y 

G/G0 

kjU)/ 

k/a)/ 

G/G0 

- -1 

k^a)/ 

k^a)/ 

_  „l/2 

Poa 

Po*1* 

q0aw 

q0a1'2 

-0.01 

3.303 

1.817 

-0.001 

1.864 

0.001 

1.365 

-0.1 

3.367 

1.835 

-0.014 

1.881 

0.010 

1.371 

-0.25 

3.480 

1.865 

-0.035 

1.908 

0.026 

1.381 

-0.5 

3.682 

1.918 

-0.072 

1.953 

0.052 

1.397 

-0.75 

3.906 

1.973 

-0.111 

1.997 

0.079 

1.411 

i 

o 

4.152 

2.032 

-0.152 

2.040 

0.107 

1.424 

-1.25 

4.422 

2.094 

-0.195 

2.082 

0.134 

1.437 

n 

4.719 

2.159 

-0.241 

2.123 

0.161 

1.448 

-2.0 

5.401 

2.299 

-0.338 

2.202 

0.216 

1.468 

-2.5 

6.215 

2.453 

-0.444 

2.277 

0.269 

1.485 

o 

CO 

1 

7.178 

2.620 

-0.559 

2.349 

0.320 

1.499 
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Table  XX.  Normalized  stress  intensity  factors  and  strain  energy  release  rates  for 
hj=a,  h^a,  v  =  0.3,  when  y  >  0. 


uniform  normal  stress 

uniform  shear 

Y 

G/G0 

kjla)/ 

k/a)/ 

G/G0 

k^a)/ 

k^a)/ 

„l/2 

Poa 

_  \n 

Poa 

q0al/2 

9oa 

0.1 

3.228 

1.797 

0.013 

1.844 

-0.010 

1.358 

0.25 

3.130 

1.769 

0.033 

1.816 

-0.025 

1.347 

0.5 

2.979 

1.725 

0.064 

1.769 

-0.050 

1.329 

0.75 

2.843 

1.684 

0.094 

1.721 

-0.073 

1.310 

1.0 

2.720 

1.645 

0.122 

1.673 

-0.095 

1.290 

H 

2.508 

1.574 

0.172 

1.577 

-0.137 

1.248 

2.0 

2.334 

1.512 

0.217 

1.482 

-0.173 

1.205 

2.5 

2.191 

1.458 

0.256 

1.389 

-0.204 

1.161 

3.0 

2.073 

1.410 

0.291 

1.300 

-0.230 

1.117 
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Chapter  5 

Conclusions  And  Future  Work 


5.1  Conclusions 

For  general  fracture  mechanics  problems,  a  half  crack-length  to  thickness  ratio 
of  approximately  unity  may  be  the  practical  limit  for  applications.  If  we  consider  an 
interface  crack  problem  that  involves  thin  films  coatings,  a  defect  in  a  solid  lubricant 
coating  as  described  in  [3]  for  example,  where  the  material  thickness  is  in  the  pm  (1.0 
x  10*6  meter)  range,  an  aspect  ratio  of  one  is  still  only  a  small  defect.  The  importance 
of  interface  cracks  in  thin  films,  therefore,  can  not  be  overstated. 

From  the  cases  of  different  geometry  and  material  constant  combinations 
computed  in  this  work,  we  may  summarize  the  following  conclusions, 

(1) .  Under  uniform  normal  stress,  the  dominant  Mode  I  stress  intensity  factor 
decreases  as  property  of  the  nonhomogeneous  material  changes  from  "soft"  to  "stiff"; 
the  nonhomogeneity  constant  changes  from  negative  to  positive.  Similarly,  the  Mode 
II  stress  intensity  factor  under  uniform  shear  shows  the  same  trend. 

(2) .  Again  speaking  of  only  the  dominant  Mode  stress  intensity  factors  under 
respective  loading,  considering  everything  else  to  be  equal,  a  decrease  in  material 
thickness,  whether  it  is  the  homogeneous  material  or  nonhomogeneous  material, 
causes  a  corresponding  increase  in  stress  intensity  factor,  we  see  from  (1)  and  (2)  two 
competing  factors  that  would  affect  stress  intensity  factors,  namely  the  degree  of 
nonhomogeneity  and  the  material  thickness. 

(3) .  Even  though  the  crack  tip  characteristics  are  of  mixed  mode  because  either 
material  constants  and  geometry  is  not  symmetric,  the  nature  of  the  loadings  is  such 
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that  the  primary  mode  is  dominant  under  respective  loading,  Mode  I  is  dominant 
under  uniform  normal  stress  and  Mode  II  is  dominant  under  uniform  shear.  The 
strain  energy  release  rate  under  respective  loadings  therefore  follows  the  trend  of  only 
the  dominant  mode  stress  intensity  factor. 

(4) .  Similar  to  their  effect  on  the  stress  intensity  factors,  the  material  thickness 
and  the  nonhomogeneity  constant  are  the  two  competing  parameters  influencing  the 
probable  angle  of  extension  of  the  interface  crack.  The  crack  extension  angle  tends 
to  turn  toward  the  material  is  "soft"  or  where  the  ligament  is  slim. 

(5) .  Effect  of  the  changes  in  Poisson’s  ratio  on  the  stress  intensity  factors  is 
small,  except  for  very  large  nonhomogeneity  constant  (y  <  -2.0).  A  Poisson’s  ratio  of 
0.3  is  good  for  practical  purposes. 

(6) .  The  most  important  information  that  can  be  obtained  from  the  crack 
opening  displacements,  if  nothing  else,  is  to  see  whether  any  compatibility  condition 
is  violated.  Because  the  crack  surfaces  are  surfaces  where  traction  boundary 
conditions  are  applied,  their  displacements  must  be  kept  in  check  to  make  sure  no 
inter-penetration  occurs.  When  it  does,  the  problem  requires  reformulation  as  a  crack- 
contact  problem. 

5.2  Remarks  on  future  research 

Computer  software  in  algebraic  manipulation  such  as  REDUCE,  MACSYMA, 
MAPLE  and  the  recently  introduced  Mathematics  which,  in  contrast  to  other  software 
and  programming  languages  that  deal  with  numerical  operations,  excels  in  doing 
symbolic  computations  and  has  broken  new  grounds  for  research  in  a  number  of 
problems  previously  considered  being  not  feasible  because  of  their  inherent  analytical 
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complexities.  Coupled  with  advances  in  computer  hardware  in  terms  of  computing 
speed  and  available  memory  as  well  as  software  development  in  user  interface,  more 
problems  are  being  attempted  with  much  efficiency.  Specific  tasks  facilitated  by  these 
symbolic  manipulators  in  helping  to  solve  research  problems  similar  to  this  study  is 
the  derivation,  simplification,  and  bookkeeping  of  very  large  number  of  terms  and 
expressions  in  the  Fredholm  kernels.  Without  such  computer  algebra  software,  it  is 
not  possible  to  solve  the  8  by  8  system  of  linear  algebraic  equations  and  to  undertake 
the  extensive  asymptotic  analysis  described  in  this  study.  Yet  what  was  utilized  in 
this  undertaking  barely  scratch  the  surface  of  the  many  faceted  features  these 
symbolic  manipulators  have  to  offer.  The  vast  potential  of  these  software  is  yet  to  be 
tapped  with  the  ever  expanding  progress  in  computing  power  that  is  reaching  new 
horizons  everyday. 

On  the  other  hand,  it  is  also  dangerous  to  assume  these  software  to  be  all 
powerful  to  solve  any  kind  of  problem  that  is  fed  into  it  without  the  user’s  careful 
planning  and  thinking  ahead.  Brute  force  is  sure  to  fail  miserably  as  problems 
attempted  get  larger  and  more  complicated,  for  sometimes  the  complexity  of  the  task 
grows  exponentially  as  the  problem  size  increases.  It  is  especially  true  when  one 
deals  with  matrix  inversion  which  requires  expanding  determinants  that  contain 
symbolic  elements.  Case  in  point  is  the  asymptotic  expansion  of  the  8  by  8  and  7  by 
7  determinants  described  in  this  work.  Without  proper  analytical  ground  work,  or 
without  correct  understanding  of  the  nature  of  the  problem,  one  has  no  choice  but  to 
rely  on  raw  computing  power  which  resulted  in  the  analysis  quickly  grinds  to  a  halt. 
Only  after  careful  planning  and  comprehension  of  the  available  computing  resources, 
as  well  as  devising  ways  to  work  around  the  limitations  could  one  secure  a 
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satisfactory  completion  of  the  task.  It  proves  once  again  that  there  is  still  no 
substitute  for  sound  analytical  reasoning.  Unlike  numerical  software,  which  always 
returns  numbers  good  or  otherwise,  symbolic  software  usually  refuses  to  work  when 
told  to  accomplish  something  that  is  beyond  its  limitations. 

In  light  of  the  capabilities  symbolic  manipulation  software  afforded  us,  a 
number  of  problems  similar  in  nature  to  what  has  been  done  in  this  dissertation 
could,  without  much  difficulty,  be  solved.  Considering  the  fact  that  the  nature  of 
materials  described  in  this  work  does  not  possess  a  distinct  interface,  but  rather  the 
material  composition  and  properties  go  though  a  transition  from  one  to  the  other. 
Also  consider  that  as  long  as  the  properties  of  the  two  materials  and  in  the  interfacial 
zone  are  continuous,  the  formulation  does  not  change  much.  An  immediate  extension 
of  the  present  problem  would  be  that  a  crack  parallel  to  the  free  surface  could  exist 
any  where  in  the  interfacial  zone. 

Cracks  of  orientations  parallel  to  the  y-axis,  with  the  crack  spanning  the 
interface,  or  the  crack  tip  terminating  at  the  interface  or  approaching  the  free  surface, 
or  an  edge  crack  are  also  prospective  problems  of  interest.  These  problems  have  the 
same  Navier’s  equations  with  constant  coefficients  as  in  this  study.  The  boundary 
conditions  will  be  somewhat  different  and  there  will  only  be  one  unknown  density 
function  to  be  solved.  Furthermore,  cracks  terminating  at  the  interface  and  crack  tip 
approaching  the  free  surface  will  have  different  characteristics  from  that  in  this  study. 
But  the  solution  of  these  problems  are  all  attainable  because  of  the  available  computer 
algebra  software  described  earlier. 

As  has  been  pointed  out  in  the  beginning  of  this  dissertation,  one  of  the 
common  causes  that  composite  materials  fail  is  the  existence  of  residual  stresses. 
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These  stresses  come  about  because  the  temperature  which  the  composite  is  processed 
differs  considerably  from  that  under  which  it  usually  operates.  This  is  also  one  of  the 
incentives  to  have  a  graded  interfacial  zone.  It  is  therefore  very  desirable  to  study  the 
fracture  mechanics  behavior  of  this  and  similar  problems  under  residual  stresses 
caused  by  thermal  mismatch.  But  until  appropriate  material  property 
characterizations  such  as  the  coefficients  of  thermal  expansion  of  the  graded  material 
are  known,  proper  assumptions  will  have  to  be  made  in  order  to  try  problems  of  this 
nature. 
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Appendix  A 
Boundary  Conditions 


In  deriving  the  system  of  singular  integral  equations  in  Chapter  2,  we  applied 
Fourier  transforms  to  Navier’s  equations  to  convert  the  coupled  PDE’s  to  ODE’s  (Eqn. 
(10)).  The  solution  to  the  ODE’s  contain  eight  undetermined  coefficients,  where  four 
are  dependent.  These  procedures  are  applied  to  both  the  homogeneous  and  the 
nonhomogeneous  materials  of  the  problem  at  hand,  which  resulted  in  a  total  of  eight 
unknown  coefficients  in  place  of  the  two  displacement  components  that  we  defined  as 
the  original  unknowns  (Eqns.  (11)  and  (18)).  A  set  of  new  unknowns  called  density 
functions  are  then  defined  as  the  x  derivative  of  the  crack  opening  displacements 
(Eqn.  (20)).  Let  the  Fourier  transforms  of  the  density  functions  be  the  new  unknown 
variables  (Eqn.  (21)).  By  taking  the  Fourier  transforms  of  the  six  homogeneous 
boundary  conditions  and  the  two  displacement  boundary  conditions  of  the  mixed 
condition,  we  can  create  a  new  set  of  eight  boundary  conditions.  This  new  set  of 
boundary  conditions,  including  six  homogeneous  condition  and  two  where  the  right- 
hand-side  are  the  new  unknown  variable  (see  Eqns.  (22)  through  (28)),  i3  what  ties  the 
eight  unknown  coefficient  functions  and  the  two  new  unknown  variables  together. 
Writing  each  of  the  new  set  of  boundary  conditions  out  explicitly,  we  obtain  an  8  by 
8  linear  system  to  solve  for  the  eight  unknown  coefficient  functions  in  terms  of  the  two 
new  unknown  variables.  The  linear  system  is  written  as  follows: 

^15^1 +^16^2 +^17^3 +®18^4  ”  (A.1) 
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**25  +  **26  ^  +  **27^3  +  ** 28  ^4  ”  0>  (A.2) 

**31^1  +  **32^2  +  **33^3  +  **34^4  "  0,  (A-3) 

a41A1+a42A2+a43A3+a44A4  «  0,  (A.4) 

**51^1  +  **52^2  +^53i<^3  +  **54^4  +  **55^-'l  +  **56^'2  +  **57^3  +  **58^-'4  ”  (A*5) 

+  ^62^2 +®63^3  +  ^64-^4  +  ^65^1  +  ^66^2 +®S7  ^3 +®68  ^4  *  (A-6) 

**71-^1  +  **72-^2  +  **73^3  +a74-^4  +  **7sCl  +^6^2  +  d77C3  +  a78C4  *  Fj,  (A-7) 

**81-^1  +  **83-^3  "^SS^l +'**86^-''2  +  **87^3  +  **88  ^4  "■^2>  (A.8) 


where  Fj  and  F2  are  the  Fourier  transforms  of  the  density  functions  that  have  been 
defined  in  Chapter  2.  The  coefficients  a  if  i,j-  1  •••■  8,  are  all  expressed  in  terms  of 
a,  k,  y,  hj  and  h2,  which  are  material  and  geometry  constants  also  defined  in  Chapter 
2.  The  actual  combining  of  terms  and  simplifying  of  expressions  with  respect  to  their 
fundamental  variable  and  material  as  well  as  geometry  constants  were  carried  out  on 
a  VAX-8300  using  the  symbolic  software  MACSYMA.  We  list  the  expressions  for  each 


element  as  follows 

a31  -  2\a\e'h'M, 

(A.9) 

a32  -  -[2A1|a|+tc-l]e-',,|o|( 

(A.  10) 

**33  -  -2  1  a ie*,la|, 

(A.11) 
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a34  -  [2/i1  ja |  -K  +  l]e*,la|, 


(A.12) 


a41  -  -2io(K-l)e'**,a,t 


(A.13) 


a42  -  ia(K-l)[JL_l+2^1]e-'*'a, 
ja 


(A.14) 


a*.  •  -2ia(K-l)«A*lal, 


(A- 15) 


-  -ia(K-l)[Sll-2h1]eA-a, 
a 


(A.  16) 


® si  ”  ^53  *  -2ia(K-l), 


(A.17) 


a52  *  1(^-1)  “ 
a 


(A.18) 


Q54  “  ~a52> 


(A.  19) 


a5i  -  -{(K  +  l)iJ3c3,9-i7(K-l)fi2e2i#-[4(K  +  5)a2+iy(K  +  l)] 


i?e‘e-4Ya2(K-ll)+f(K  +  l)} _ ilK~1) _ , 

8a[i2e‘8-”y(2  -k)] 


(A.  20) 


ase  ■  -((K  +  l)fl3c'3ie-t7(K-UR2e'2,e-[4(K  +  5)a2+i72(K  +  l)] 


Re  "ie -4ya2(K  - 11)  +t3(k  + 1)} 


U  K-l) 

8a[Re~‘6  -y{2  -k)]’ 


(A.21) 
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(A.22) 

(A.  23) 

(A- 24) 

(A- 25) 

(A.26) 

(A.27) 

(A.  28) 

(A.29) 

(A.30) 

(A.31) 
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(A.32) 


(A.33) 

(A.34) 

(A.35) 

(A.36) 

(A-37) 

(A.38) 

(A.39) 

(A.40) 

(A-41) 

(A.42) 
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-(k-  l)i?2e2ie  +  4(K  +  l)a2+y2(K- 1) 
4[i?e‘e  -y(2  -k)] 


(A.  43) 


«85  *  a86  "  °87  *  a88  *  * «.  (A.43) 

where  R  and  9,  both  defined  in  Chapter  2,  are  shown  again  as  follows, 

R  » 

(A.49) 
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Appendix  B 
Asymptotic  Expansions 

B.l  Introduction 

In  the  derivation  of  the  system  of  singular  integral  equations,  an  important 
step  is  finding  the  asymptotic  values  of  the  integrands  of  some  infinite  integrals 
(Eqns.  (47)  to  (50)).  The  asymptotic  analysis  is  necessary  for  the  following  two 
reasons, 

1.  The  leading  term  of  the  integrands  in  asymptotic  expansion  gives  the 
Cauchy  kernel  in  the  singular  integral  equations. 

2.  Subsequent  terms  in  the  expansion  facilitate  computational  efficiency 
when  we  numerically  solve  the  SEE’s.  The  more  terms  that  can  be 
asymptotically  extracted,  the  less  effort  is  required  for  the  numerical 
computation. 

The  major  task  of  the  analysis  is  the  asymptotic  expansion  of  the  8  by  8 
determinant  and  eight  7  by  7  cofactors  which  form  the  integrands  in  the  infinite 
integrals.  These  cofactors  and  the  determinant  constitute  the  solution  of  the  8  by  8 
linear  system  as  depicted  by  Eqn.  (29).  The  leading  terms  of  i  =  1  •—  A,j  =  1, 

2  in  Eqn.  (30)  are  desired. 

B.2  Algebraic  vs.  numeric  operations 

One  may  wonder  why  in  solving  the  linear  system  (29),  the  inefficient  Cramer’s 
rule  was  used  in  lieu  of  more  efficient  methods  (see  Eqn.  (30)  and  (31)).  The  answer 
to  this  has  to  do  with  the  inherently  different  nature  of  algebraic  manipulations  and 


numeric  operations.  Consider  the  problem  of  inverting  an  n  by  n  matrix.  What  we 
learned  from  numerical  analysis  tells  us  that  the  problem  requires  operations  on  the 
order  of  n3  [31].  But  this  kind  of  analysis  does  not  apply  when  we  are  dealing 
with  algebraic  entities.  Consider  now  the  number  of  terms  in  the  determinant  of  a 
fully  occupied  n  by  n  matrix  whose  elements  are  the  symbols  atJ.  The  determinant, 
when  fully  expanded,  which  is  necessary  in  computing  the  inverse,  has  n!  terms.  For 
an  8  by  8  matrix  with  each  element  being  a  single  symbol,  the  fully  expanded 
determinant  has  40,320  terms,  and  this  is  vastly  different  from  the  number  of 
operations  (8s  =  512)  needed  to  invert  the  matrix  numerically.  For  the  problem  at 
hand  where  each  element  in  the  8  by  8  matrix  is  a  complex  expression  rather  than  a 
single  symbol,  the  number  of  terms  in  the  expanded  determinant  is  even  more.  The 
complexity  of  determinant  expansion  will  quickly  make  the  asymptotic  analysis  come 
to  a  halt  if  brute  force  determinant  expansion  was  used. 

A  straight  forward  approach  would  be  to  asymptotically  expand  each  element 
in  the  matrix  to  compute  the  determinant  and  cofactors  by  direct  expansion.  One  soon 
finds  that  such  a  procedure  breaks  down  quickly  as  matrix  size  and  number  of  terms 
expanded  in  each  element  increases.  The  reason  for  the  break-down  is  obvious.  As 
the  matrix  size  increases,  even  if  there  is  only  one  term  in  each  element,  the  number 
of  terms  in  the  expanded  determinant  increases  exponentially  as  explained  earlier. 
When  there  are  more  terms  in  each  element,  and  each  term  is  in  turn  complex 
expressions,  it  soon  becomes  more  than  even  a  mainframe  computer  can  handle.  On 
a  VAX-8530  in  which  20,000  pages  of  dynamic  memory  was  allocated,  the  computation 
breaks  down  with  a  two  term  expansion  in  each  element  of  an  8  by  8  matrix  to 
analytically  expand  its  determinant. 
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The  successful  asymptotic  expansion  of  the  8  by  8  determinant  and  its  cofactors 

to  as  many  terms  as  the  aforementioned  restraint  would  allow,  lies  in  recognizing  that 

out  of  the  nominally  required  astronomical  number  of  multiplications,  only  a  relevant 

few  need  to  be  carried  out  to  produce  the  leading  terms  in  the  determinant.  Those 

expansions  not  carried  out  contribute  only  to  low  order  terms  which  can  be  ignored. 

In  what  follows  each  relevant  element  of  the  8  by  8  matrix  undergoes  a  12-term 

asymptotic  expansion.  These  expansions  in  turn  produce  a  12-term  asymptotic 

D 

polynomial  out  of  (see  Eqns.  (51)  through  (54))  which  gives  a  11 -term 

asymptotic  analysis  in  the  integrands  of  infinite  integrals  (Eqns.  (47)  through  (50)). 
The  memory  constraint  is  the  only  factor  that  controls  the  number  of  terms  that  can 
be  asymptotically  expanded  in  each  element,  therefore  it  also  determines  the  number 
of  terms  that  can  be  expanded  in  the  integrands  of  Fredholm  kernels. 


B.3  Asymptotic  expansion 

To  asymptotically  expand  the  8  by  8  determinant  and  eight  7  by  7  cofactors, 
we  first  examine  the  asymptotic  behavior  of  the  elements  that  make  up  these 
determinants.  Referring  to  expressions  of  each  element  in  Appendix  A,  we  note  that. 


R 


N 


-  2a  +  0(0), 

1  +  K 


0  -  Itan-H-iST-  , 

2  y2+4a2‘v] 


ill  >-Iv 

1+K  2  \1 


~  -  +  o(JL), 

1+K  a  a2 


(B.l) 


a  — » 

therefore, 
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Rei6  -  2  a  +  0(0), 

(B  St) 

Re'1*  -  2a  +  0(0),  as  a  -»  <». 

As  have  been  defined  in  Chapter  2,  O(-L)  indicates  polynomials  with  leading  term 

a2 

of  degree  -2  in  a,  and  so  on.  With  these  observations,  we  found  that  a1S)  aJ6 ,  and  a25, 
a26  contain  a  term  like  eaK>  in  their  asymptotic  behavior.  Similarly,  elements  a33> 
a34,  and  a43  ,  au  have  a  term  like  e°*‘  .  We  bear  in  mind  that  in  expanding  a 
determinant,  no  two  elements  that  belong  to  the  same  column  or  row  may  appear  in 
the  form  of  product  in  any  expansion.  It  is  apparent  that  the  dominant  term  in  both 
the  numerator  and  denominator  must  carry  a  term  like  e  2a(h'  * h%>  so  that  the  leading 
terms  are  derived  from  their  quotient.  Since  the  denominator,  the  determinant  of  the 
8  by  8  matrix,  will  always  carry  the  term  e2a(hx  * A>)  ,  any  term  in  the  numerator 
without  e2a<h'  *  will  only  contribute  to  lower  order  terms  in  the  rational 
expression.  Our  objective  then,  becomes  to  identify  all  relevant  elements  that 
contribute  to  this  dominant  term  in  each  8  by  8  or  7  by  7  determinant  expansions  and 
ignore  the  rest  as  they  are  insignificant  to  the  analysis. 

The  foregoing  analysis  is  the  key  to  the  successful  asymptotic  analysis  of  the 
problem.  It  reduces  what  appeared  to  be  an  insurmountably  complicated  task  to  a 
manageable  job  that  can  be  better  explained  by  rewriting  the  8  by  8  matrix  in  a 
different  way  as  shown  in  Eqn.  (B.3).  Following  the  rule  in  expanding  a  determinant 
as  stated  earlier,  we  observe  that  the  dominant  term  in  the  determinant  must  include 
two  elements  from  columns  3,  4  and  two  from  columns  5,  6  with  the  exponential  term. 
Another  way  of  getting  the  dominant  terms  is  through  row-wise  expansion,  in  which 
case  the  leading  terms  of  the  determinant  must  include  two  elements  from  rows  1,  2 
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0  0  0  0  -«„«“**  7  7 

0  0  0  0  -asseah'  - aseeah *  7  7 

I  I  0( a)eah'  0{ct)eah'  0  0  0  0 

7  7  0{a.)eah'  CKoOe^1  0  0  0  0 

(B.3) 

®51  a52  I  I  I  I  a57  a  58 

a6i  as2  7  7  7  7  a67  agg 

®71  ®72  I  I  I  I  G77  fl7g 

a8i  0  7  0  7  7  a87  aM 

and  two  from  rows  3,  4  that  have  the  exponential  term.  With  these  elements  in  the 
product,  the  rest  will  have  to  come  from  those  elements  explicitly  written  in  Eqn. 
(B.3),  observing  the  rule  that  only  one  element  in  a  column  or  row  may  participate  in 
forming  the  product.  The  remaining  elements  are  either  zero  occupied  or  do  not 
contribute  to  the  leading  terms  of  our  analysis  and  are  marked  as  "7". 

Having  established  the  ground  rule  of  obtaining  the  leading  terms  of 
asymptotic  analysis  in  the  determinants,  we  can  now  explicitly  write  the  exact 
expressions  to  get  the  leading  terms  of  the  8  by  8  determinant  and  7  by  7  cofactors. 
Note  that  they  are  written  in  terms  of  truncated  asymptotically  expanded  expressions 
of  only  those  elements  associated  with  products  that  give  the  dominant  terms.  Let  us 
denote  the  truncated  asymptotic  expressions  of  elements  as  a\j,  and  those  of  the 
determinant  and  cofactors  as  d  and  di}  respectively.  We  have 
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(B.4) 


rfjj  ■  (Ogs'Ogg “Ogj'Ojg)  *(034*043  ”033*044)  *[052* 


'a67-a88  ”Og7*Ogg)  ”0g2*(058*0g7  ”057"Ogg)], 


,  1  /  /  '  \  ,  I  i  /  /  ,  r  /  ,  t  /  /  /v 

”'055*flgg  “Ogg'fljg)  *(034*043  ”033*044)  *L052\Qg7*Q7g  ”077"Ogg)  ' 


l,l  /  /  /  v  /  ,  /  /  /  /  n 

®57  '^62*^78  ~®72’a68'  +®58  va62’®77  ”072"Og7)J, 


<f2l  *  -(o5s*Ogg  “Ogj-fljg)  *(034*043 -033*044)  *[a5i^ag7‘flg8”fl87'Og8)' 


1,1  /  /  /  .  /  ,  /  /  /  /  ,i 

fl57*'®61*®88  "Ogi’flgg)  +  Ogg  *(Ogi‘Og7  ”08i'Og7)j, 


.  ,  /  /  /  /  .  ,  /  /  /  >  \  r  >  j  '  '  / 

022  *  va55*®6« -a«8*056' *'®34‘®43 -033*044) 'Lfl5i'(fl67’a78  ~®77‘a68' ' 


1,1  /  /  /v  /  ,  /  /  /  /  ., 

057*(fl6i*fl7g  On’aeg)  +o58  *(asl*aT!  o71*a67)J, 


c?  •  ”(oS5*Ogg  o68*ass) -(034*043  O33 -044)* {o8j- [Og2‘(Og7’a7g ~ o77*agg)  + 


/  ,  I  /  _/  /  \-,  /  r  /  /  /  /  /  />.  /  ,  /  / 

072*'0g7*0g8  —  Og7*Ogg/J  -a52ia61 -(087*078  -077-088)  +a 72‘(08i‘Q78- 


(B.5) 


(B.6) 


(B.7) 


t  /  v  /  /  /  I  /  /  /  r  /  ,  /  /  /  /  v 

07i*088)  +  Ogg *( Ogj *077  _fl7i*fl87)J  +  057*l”Og2*(fl81*078  ”071*088)  ' 


(B.8) 


/  /  r  /  /  /  /  /  /\  /  /  /  / 

072*(Ogi  *Ogg  -Ogi’flgg))  —Ogg'L” flg2<(Ogj*fl77  ” fl7l’Og7) -a72*(ag1-a87 ' 


Ogi*Og7)]} , 


It  is  amazing  that  there  are  only  54  terms  required  to  obtain  the  leading  terms  of  the 
8  by  8  determinant,  out  of  the  possible  40320  (8!),  as  shown  in  Eqn.  (B.8).  Of  course 
the  number  of  terms  will  be  something  less  if  full  determinant  expansion  was  carried 
out  because  there  are  18  zero  occupied  elements  in  our  8  by  8  determinant  (see  Eqn. 
(B.3)).  But  the  reduction  in  the  number  of  terms  needed  to  compute  the  leading  terms 
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of  the  determinant  paves  the  way  for  a  promising  prospect  that  perhaps  more  terms 
can  be  expanded  from  our  asymptotic  analysis.  It  is  also  worth  noting  that 
expressions  for  d3I,  d32,  d41,  and  d42  need  not  be  written.  The  reason  for  this  is 
apparent  when  we  notice  that  the  corresponding  cofactors  D31,  D32>  D41,  and  D42  were 
obtained  from  the  8  by  8  matrix  by  deleting  column  3  or  4.  In  any  case  one  of  the 
dominant  terms  eah'  was  lost  thus  making  these  cofactors  of  much  lower  order  than 
that  of  determinant  D.  Therefore  the  resulting  rational  polynomials  are  analytic  at 
a  =  oo  when  such  lower  order  cofactors  appear  in  the  numerators.  In  other  words, 

— — ,  P.—,  and  P42  do  not  contribute  to  Cauchy  kernel  or  any  of  the 

D  D  D  D 

leading  terms  in  the  asymptotic  analysis  that  aids  numerical  efficiency,  but  become 
part  of  the  massive  Fredholm  kernels. 

Due  to  the  constraint  in  memory  allocated  on  the  VAX-8530,  there  are  only  so 
many  leading  terms  we  can  expand  out  of  the  asymptotic  analysis.  We  shall  now  try 
to  establish  the  relationship  between  the  number  of  leading  terms  in  the  overall 
asymptotic  analysis  with  that  of  each  element  in  the  matrix.  Let  -n  be  -  lowest 
power  that  can  be  asymptotically  expanded  in  a  particular  element.  We  denote  the 
range  of  terms  that  can  be  expanded  within  the  constraint  for  an  element  as 
0(<x,m,-n),  where  m  is  the  power  of  the  leading  term  and  a  is  the  variable  of  the 
asymptotic  expansion.  We  can  write  number  of  terms  expanded  for  each  significant 
element  and  exact  asymptotic  expression  for  some  elements  in  Eqn.  (B.3)  as  follows, 


a™-  0(a,l)  -  -2  a, 


(B.9) 


a'34  -  O(a,l,0)  ■  2/ixa-K  +  l, 


(B.10) 
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a  «  “  0( a,l)  "  -2ia(K-l), 


a 'u  -  O(a,l,0)  -  i(K-lX2A1a-K  +  l), 


a  j!  -  0(a,l)  •  -2ia(»c-l), 


a  52  -  O(a,0)  «  ^k2-!), 


a'61  -  0(a,l)  -  2  a, 


a  62  »  O(a,0)  -  1-k, 


a  71  -  0(a,l)  -  -a, 


*  a87  -  a'se  -  0(a,l)  -  ia, 


_  /  _  1  _1  .  .  „  .  .  .  C65,-l  .  C58,-2  ,  C55,-3  .  C55,-4  .  ...  .  C55,-» 

a  55  *  0(a,l,-n)  -  c551a+c550  + - +  — +  — _  +  — _  +  •••  +  — — , 

a  a*  or  a4  a" 


_/  .  nin  1  _  „  „  ,  -  ,  C86,-l  C56,-2  .  C56,-3  .  C5€.-4  .  ...  ,  C56,-« 

a  56  -  0(a,l,-n)  -  c56,a+c56n  + - + — ^-  +  — =-  + — — +  -  +  — — , 


a  a2  a3  a4 


a" 


„/  -  ru~  1  _  ,  „  .  „  .  C57,-l  C57,-2  ,  CJ7,-3  C57,-4  .  ..  ,  C57,-n 

a  57  -  U(a,l,-n)  -  cs71a+c570  + - +  — •  + — :_  + — i_  +  -+ — l_, 

a  a2  a3  a4  a" 


_/  .  nin  1  ~  „  .  C58,-l  .  C58,-2  ,  C58,-3  .  C58,-4  .  ...  .  C58,-» 

a  58  C/(a,l,  Tl)  c58,l®  +C58.0  + - + - TT~  + - 7~  + - 7“  +  + - * 


a  a2  a3  a4 


a” 


_/  _  /')/_  1  ^  „  .  C65,-l  ^  C65,-2  .  C65,-3  C65,-4  .  .  ,  C65,-i> 

a  6S  -  U(a,l,-n)  -  c68ila+c660  + - +  — l_  + — —  + — —  +  “•+ — — , 

a  a2  a3  a4  a" 


a72  -  O(a,0)  -  k, 


(B.ll) 

(B.12) 

(B.13) 

(B.14) 

(B.15) 

(B.16) 

(B.17) 

(B.18) 

(B.19) 

(B.20) 

(B.21) 

(B.22) 

(B.23) 

(B.24) 
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where  cij  k  represents  the  coefficient  of  the  ot*  term  of  the  truncated  polynomial  a\j. 
Note  that  elements  a’33>  a’34>  a’^  and  a’^  were  written  without  the  term  e0*1  .  Also 
a’ is  and  a’ie  differ  from  a’65,  and  a’^  by  the  factor  -eah>  ,  as  can  be  seen  from  Eqns. 
(A.32),  (A.33)  and  (B.2),  the  latter  is  written  instead.  For  the  same  reason  a’S5,  and 
a’56  are  written  instead  of  a’2S,  and  a’26.  The  reason  that  both  eahl  and  eaAj  are 
excluded  from  the  analysis  is  that  they  factor  out  since  they  appear  in  both  the 
numerator  and  the  denominator  when  RRl  is  divided.  Eqns.  (B.9)  through  (B.18) 
are  those  elements  whose  exact  representations  are  simple  polynomials.  Therefore, 
no  asymptotic  expansions  are  needed.  They  art  lirectly  constructed  from  Appendix 
A.  Eqns.  (B.19)  through  (B.29),  however,  are  those  elements  that  have  complex 
rational  polynomial  expressions.  We  shall  apply  asymptotic  expansion  to  these 
elements,  and  their  exact  expressions  up  to  the  lowest  power  expandable  will  be 
shown  later  in  the  appendix. 
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If  one  tries  to  substitute  Eqns.  (B.9)  through  (B.29)  into  Eqns  (B.4)  through 
(B.8)  now,  considering  that  the  necessary  asymptotic  analysis  has  been  done  on  the 
elements  shown  in  Eqns.  (B.19)  through  (B.29)  and  all  the  coefficients  ctJ  k  are  known, 
one  still  can  not  expect  to  extract  the  most  leading  terms  out  of  the  determinant  and 
cofactors.  The  reason  is  as  follows.  Remember  that  a’tJ  are  truncated  polynomials, 
therefore  any  operation  on  them,  especially  multiplications  and  divisions,  will  produce 
nonsignificant  terms.  If  steps  are  not  taken  to  keep  only  the  significant  terms  and  to 
throw  away  the  excess  baggage,  before  very  long  these  nonsignificant  terms  will 
accumulate  to  a  monstrous  proportion  and  defeat  the  economizing  that  we  have 
accomplished  earlier.  Therefore,  we  now  establish  the  following  rules  as  regard  to 
operations  on  these  truncated  polynomials 

0(a,nlf  -n2)+0(a,n3,  -n4)  -  0(a,max[n1,n3l,maxl-n2, -nj), 


0{a,nv-n2)-0(a,n3,-n4)  -  0(a,max[n1,n3],max[-n2, -nj), 


0( a,  nv  -n2)  -0(a,  n3,  -n4) 


0{  a,  +  n3,max[n3  -  n2,  -  nj), 


0(a,nv  -n2) 
0(a,  n3,  ~fi^) 


0{a,n1  -n3,  -n2-n3). 


(B.30) 


We  can  now  substitute  all  the  expressions  in  Eqns.  (B.18)  through  (B.29)  into 
Eqns.  (B.4)  through  (B.8)  to  obtain  the  significant  expanded  terms  in  the  cofactors  and 
determinant.  The  resulting  range  of  significant  terms  are 
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dn  -  0(a,6,4-n), 


d12  ■  0(a,6,4-n), 

du  -  0(a,7,5-n),  <B*31> 

d22  -  0(a,7,5-n), 

-  0(a,7,5-n), 

The  actual  asymptotic  expansion  was  done  on  VAX-8530  using  the  symbolic 
software  MAPLE.  It  gives  n  =  12,  which  is  the  limit  to  which  the  asymptotic 
expansion  of  those  elements  in  Eqns.  (B.19)  through  (B.29)  can  be  carried  out. 
During  the  asymptotic  expansion  process  for  those  elements  in  Eqns.  (B.19)  through 
(B.29),  it  was  found  that 


*  C55,l0t~C55,0  + 


*  C56,ia  C56,0  + 


"  “C65,ia  +C66,0 


*  C66,ia  +C66,0 


CS5.-1  _ 

CSS,-2 

C56,-3  _ 

C®5.-4  4. 

(B.32) 

a 

a2 

a3 

a4 

a" 

C  56,-1  _ 

C56,-2 

C56,-3  _ 

C*6.-4 

(B.33) 

a 

a3 

a4 

a" 

_  C65,  -1 

j.  C65,-2 

_  C65,-3 

4.  C65.-*  . 

—  -*-(— 1)B  Cs5‘'*, 

(B.34) 

a 

a2 

a3 

CL* 

a" 

_  C 66,-1 

x  C66,-2 

_  C66,-3 

+  C66,-4  _ 

....+(_!)»  c«*-\ 

(B.35) 

a 

a2 

a3 

a4 

a" 

As  a  result,  we  need  to  show  only  the  expanded  terms  of  elements  ass,  a6S>  a77, 
a78.  They  are  as  follows, 
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Re[a's6\  «(— (k  +  1)  -  -I1k(k  +  1)  +  -J— -(k3  -3k -2)  -  t-.,  (k^-k3  - 

2  4a  8a2  16a3 

5k2-k  +  4)  +  _J1_(k5-2k4-6k3+4k2  +  13k  +  6)  -  - t - 

32a4  64a5(K  +  l) 

(k7  -  2k6  -  9k*  +  6k4  +  33k3  +  18k2  -  2  Ik  -  32)  +  _L_(k7-4k6- 

128a6 

5k5  +  22k4  +  23k3  - 32k2  - 51k - 18)  - _ t _ (k1 °-3 k9- 

256a7(K  + 1)2 


12k8  +  22k7  +  78k6  -  18k5  -  222k4  - 198**  +  63k2  +  225k  +  252)  + 


_ t _ (k9  -  6k8  +  44k6  -  6k5  -  144k4  -  56k3  +  180k2  +  189k  +  54)  - 

512a8 

_ tl _ (k13  -  4k12  -  14k11  +  46k10  +  125k9  -  170k8  -  660k7  - 

1024a9(K  +  l)3 

(B.36) 

120k6  +  1425k8  +  1820k4  +  262K3  -  1380k2  -  1709k  -2192)  + 


_ tl — (k1’  -  8k10  +  9k9  +  62k8  -  94k7  -  264k6  +  250k5  +  724k4  -  3k3  - 

2048a10 

864k2  -675k  -162)  -  _ tl _ (k16- 5k15 -15k14 +  77^®  + 

4096a11(K  + 1)4 


161k12  -483k11  -  1253k10  +  1015k9  +  5355k8  +  2975k7  -8323k6  - 


15015k5  -  7175k4  +  6347k3  +  13341k2  +  11437k  +  21052)  + 


_ tl _ (k13  -  10k12  +  22k11  +  68k10  -  245k9  -  262k8  +  1060k7  + 

8192a12 


1016k6  -  2201k5  -  3030k4  +  1062k3  +  3780k2  +  2349k +486)} 
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Imia'A  -  2a  -  1<k  +  1)  +  ^(k2 - k - 2)  -  jLtic3- 2k2 -3k  +  2)  + 
2  4a  8a2 

— ^—(^-Sk3 -S^  +  Tk  +  S)  +  - — - (k6-3k5-6k4  + 

16a3  32a4(K  +  l) 

12k3+21k2-3k-20)  +  —^—(k6 -51c5  +  22K3  +  K2 -33k- 

64  a6 

18)  -  _ t _ (k9  -  4k3  -8k7  +  30k6  +  48K6  -  66k4  - 

128a6(K  +  l)2 


156k3  -42k2  +  103k +  154)  +  — £  .(k8  -7k7  +  7K6  +  37K5  - 

256  a7 


43k4  - 101**  +  45k2  +  135k  +  54)  -  _ £ _ (k12-5k“- 

512  a8(K  +  l)3 


9k10  ♦  55k9  +  70k8  -  240k7  -  420k6  +  300k6  + 1  125k4  +  695k3  - 

435k2  -  885k  - 1270)  +  ^--.(k10  -  9k9  +  18k8  +  44k7  -  (B.37) 

1024a9 

138k6 -  126k5  +  376k4  +348K3  -351k2  - 513k  - 162)  - 


_ I _ (k16  -  6k14  -  9k13  +  86k12  +  75k11  -  558k10  - 

2048a10(K  +  l)4 

695k9  +  1710k8  +  3645k7  -  670k6  -  7653k6  -  7362k4  + 

185k3  +  6258k2  +  5823k  + 1 1638)  +  — £! — (k^-IIk11* 

4096  a11 

33k10 + 35k9  -  280k8  +  18k7  +  1042k6  -26k6  -  2175k4  - 


855k3  +  1917k2  +  1863k  +  486)  -  - £! - (k18- 

8192a12(K  + 1)6 

7k17  -  8k16  +  122k16  +  52k14  -  1036k13  -  770k12  +  4970k11  + 


7084k10  -  10780k9  -  30338k8  -  7294k7  t  47824k6  +  65744k6  + 


22052k4  -  32186k3  -  56905k2  -  29345k  - 116376), 
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Re[a' 56]  -  -Re[a'bb], 


(B.38) 


/mta'jj  -  lm[a'bb\ 


Re[a'bb\  -  -2a  +  -I( k-1)  -  - t _ (k3-2k2-3k  +  4)  +  -^L-Ck3 -4k2 + 

2  4a(K  +  1)  8a2 

k  +  6)  -  _ t _ (k6  -  3k5  -  6k4  +  12k3  +  21k2  -  k  -  28)  +  t 

16a3(K  +  l)2  32a4 

^-61^  +  6^  +  16k2 -15k -13)  -  _ t _ (k9-4k8-8k7  + 

64a5(K  +  l)3 

30k6  +  48k5  -  66k4  -  156k3  -40k2  +  79k +  212)  +  . J__ (k7-8k6  + 

128  a6 

15k5  +  22k4  -  65k3  -  36k2  +  81k  +  54)  -  _ t _ (k^-Sk11- 

256a7(K  + 1)4 

9k10  +  55k9  +  70k8  -240k7  -420^  +  300^  +  1125*4  +  697**  -483k2  - 

591k -1784)  +  — t — (k9  -  10k8  +  28k7  +  16k6  -  154k5  +  28k4  + 

512  a8 

348k3 - 351k - 162)  -  _ tl - (k16- 6k14 -9k13  +  86^*  + 

1024a9(K  +  l)5 

75k11  -  558k10  -  695k9  +  1710k8  +  3645k7  -  670k6  -  7653k5  -  7360*4  +  (B.39) 

105k3  +  7  150k2  +  2  135k)  +  tl — (k11- 12k10 +  45^-10^- 

2048a10 

270k7  +  288k6  +  754k5  -  780k4  -  1395k3  +  540k2  +  1377k  +  486)  - 

_ tl _ (k18  -  7k17  -  8k16  +  122k15  +  52k14  -  1036k13  - 

4096a11(K  +  l)6 

770k12  +  4970k11  +  7084k10  -  10780x^-30338^  -  7294k7  + 


47824k6  +  65746k5  +  21932k4  -  30086k3  -  71965k2  +  17725k  - 

169596)  +  — tl — (k13  -  14k12  +  66k11  -  64k10  -  385k9  +  858k8  + 
8192a12 

988k7  -3156k6  -2097k5  +  5670k4  +4482^-3888^  -  5103k  - 


1458), 
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Im[a'6S]  -  (-X(k-I)  -  Sfl—(y? -k-2)  +  - t - (k^-x^-S^-k  +  S)  - 

2  4a  8a2(K  + 1) 

_ I?_(k4-3k3-3k2  +  7k  +  6)  + - t - (k7  -2  k6 -9k5  + 

16a3  32a4(K  +  l)2 

6k4  +  33k3  +  18k2  -  19k - 44)  -  ^-(k^-Sk5  +  22K3  +  K2  - 

64  a5 

33k -18)  +  _ t _ (k10-3k9-12k8+22k7  +  78k6- 

128a6(K  +  l)3 

18k5  -  222k4  -  198k3  +  65k2  +  193k  +350)  -  _  J—  ;  k8  -  7k7  + 

256a 

7k6  +  37k5-43k4-1011£‘  +  45k2  +  135k +  54)  -  _ £. _ 

512a8(K  +  l)4 

(k13  -  4k12  -  14k11  +  46k10  +  125k9  -  170k8  -  660k7  -  120k6  + 

1425k5  +  1820k4  +  264k3  -  1440k2  -  1263k  -3114)  -  t°  — 

1024a9 

(k10  -  9k9  +  18k8  +  44k7  -  138k6  -  126k5  +  376k4  +  348k3  -  351k2  - 


513k -162)  +  - tl - (k16  -  5k15  -  15k14  +  77k13  + 

2048a10(K  +  l)5 

161k12  -  483k11  -  1253k10  +  1015k9  +  5355k8 + 2975k7  - 


8323k6  -  15015k8  -  7173^+6251^  +  14601k2  +  5413k  + 

30466)  -  _ tl _ (k12  -  11k11  +  33k10  +  35k9  -  280k8  +  18k7  + 

4096a11 

1042k6  -26k5  -  2175k4  -  855k3  +  1917k2  +  1863k  +486)  + 

_ tl _ (k19  -  6k18  -  15k17  +  114k16  +  174k15  -  984k14  - 

8192a12(K  +  l)6 

1806k13  +  4200k12  +  12054k11  -3696k10  -41118k9  -37632k8  + 


40530k7  + 1  13568k6  +  87800k5  -  10414k4  -  83451k3  ■ 


132210k2  -  15619k  -  319176)) 


3  -k 

^  1  +K 


(B.40) 
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Re[a  W  - 

lm[a‘^\  -  -/mta'gj], 


(B.41) 


Beta']  -  -a  -  1{  k-2)  -  _ £ - (k3 -2k2 -3k  +  2)  -  -lltK3-^* 

2  4a(K  +  1)  8a2 

k  +  6)  -  _ _ (k6  - 3k5  -  6k4  +  12k3  +  21k2  -  3k  -  20)  -  _J__ 

16a3(K  + 1)2  32  a4 

(k5 - 6k4  +  6K3  +  16k2 -  15k - 18)  -  _ £ _ (k9-4k8-8k7  + 

64a5(K  + 1)3 

30k6  +  48k5  -  66k4  -  156k3  -  42k2  +  103k  + 154)  -  -.Y— (k7  -8k6  + 

128a6 

15k5  +  22k4  -  65k3  -  36k2  +  81k  +  54)  -  _ £ - (k^-Sk11- 

256a7(K  + 1)4 

9k10  +  55k9  +  70k6  -  240k7  -  420k6  +  300k5  + 1  125k4  +  695*  -  435k2  - 

885k  - 1270)  +  — £ — (k9  -  10k8  +  28k7  +  16k6  -  154k5  +  28k4  + 

512a8 

348k3 -351k -162)  -  - H - (k15- 6k14 -9k13  +86^*  + 

1024a9(K  + 1)5 

75k11  -  558k10  -  695k9  +  1710k8  + 3645k7  -  670k6  -  7653k6  - 

(B.42) 

7362k4  +  185k3  +  6258x^+5823^  -  £1 (k11  -12k10  +  45k9- 

2048  a10 

10k8  -  270k7  +  288k6  +  754k5  -  780k4  -  1395k3  +  540k2  +  1377k  + 

486)  -  _ £1 _ (k18 -7k17 -8k16  +  122k15  +  52k14- 

4096au(K  +  l)6 

1036k13  -  770k12  +  4970k11  +  7084k10  -  10780k9  -  30338k8  - 

7294k7  +  47824k6  +  65744k5  +  22052k4  -  32186k3  -  56905k2  - 

29345k)  -  — £1 — (k13  -  14k12  +  66k11  -  64k10  -  385k9  +  858k8  + 

8192  a12 

988k7  -  3  156k6  -  2097k5  +  5670k4  +  4482k3  -  3888k2  -  5103k  - 
1458), 
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Im[a  78] 


-  {-Ik  -  J^.(k2-k-2)  -  _ £ _ (k4  -k3  -5k2-k  +  4)  - 

2  4a  8a2(K  + 1) 

.  ^(1^-3 tfSyf  +  lK  +  G)  -  _ £ _ (k7  -2k6-9k®  + 

16a3  32a4(K  +  l)2 

6k4  +  33k3  +  18k2  -  21k  -  32)  -  _7L(k6-5k5+22k3+k2- 

64  a5 

33k -18)  -  _ £ - (k10  -  3k®  -  12k8  +  22k7  +  78k6  - 

128a6(K  +  l)3 

18k5  -  222k4  -  198k3  +  63k2  +  225k  +  252)  -  -  (k8  -  7k7  + 

256  a7 

7k6  +  37k5  -  43k4  - 101*’  +  45k2  +  135k  +  54)  -  J_ 

512a8(K  +  l)4 

(k13  -4k12  -  14k11  +  46k10  +  125k9  -  170k8  -  660k7  -  120k6  + 

1425k5  +  1820k4  +  262k3  -  1380k2  -  1709k  -2192)  -  £. ° — 

1024a9 

(k10  -  9k9  +  18k8  +  44k7  -  138k6  -  126k5  +  376k4  +  348k3  -  351k2  - 


513k -162) 


Y11 


-(k16  -  5k15  -  15k14  +  77k13  + 


(B.43) 


2048alo(K  +  l)5 
161k12  -  483k11  -  1253k10  +  1015k9  +  5355k8  +  2975k7  - 


8323k6 - 15015k® -7175k4  +  6347k3  +  13341k2  +  11437k  + 

21052)  -  — £1 _ (k12  -llKn  +  33k10  +  35^-280^  +  18k7  + 

4096a11 

1042k6  -26k5  -  2175k4  -  855k3  +  1917k2  +  1863k  +486)  - 

- £1 _ (k19  -  6k18  -  15k17  +  114k16  +  174k15  -  984k14  - 

8192a12(K  +  l)6 

1806k13  +  4200k12  +  12054k11  -  3696k10  -  41118k9  -  37632k8  + 


40530k7  + 1 13568k6  +  87798k5 -  10274k4 - 86271k3 - 


109230k2  -  65051k  -  2 17776)) 


N 


3  -k 

1  +  K 
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Reia'v]  -  Re[a'77], 

(B.44) 

lm[a'a7]  •  -lm[a'77]. 

While  substituting  Eqns.  (B.9)  through  (B.29)  into  Eqns.  (B.4)  through  (B.8), 
some  of  the  coefficients  in  the  leading  terms  of  expansion  become  zero  when  the  actual 
computation  was  done.  As  a  result,  instead  of  seeing  the  power  of  leading  terms  as 
shown  in  the  second  parameters  in  Eqn.  (B.31),  we  obtain  the  following 

dn  -  0(a,3,-8), 

d12  •  0( a, 3, -8), 

d21  -  CXa.4,-7),  ®-*5) 

d22  •  0(a,4,-7), 
d  -  0( a, 4, -7). 

Using  the  rules  set  forth  in  Eqn.  (B.30),  the  final  asymptotic  expansions  are  as 
follows, 

~  -0(a, -1,-12), 

d 

-0(«, -1,-12), 

a 

(B.46) 

-Hi  -  O(a,0,-ll), 

a 

4?  *  0(a, 0,-11). 
a 

Written  explicitly,  they  are 
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-  - _ — (k - 1)  + _ 1 _ (k-3)- - 1 - (3k*-k-6)  + 

d  a(x  +  l)  4a2(K  +  l)2  8a3(K  +  l)3 


_ I _ (k3  -  8k2  +  17k  + 10)  - _ I _ (3k4  -  21k3  +  17k2  +  5k  • 

16a4(K  +  l)4  32a5(K  +  1)5 


4)  + _ £ _ (k5  -  17k4  +  76k3  -  72k2  -  173k  -  71)  -  —£ _ ■ 

64a6(K  +  l)6  128a7(K  +  l)7 


(3k6  -  50k5  +  159k4  -4k3  -  127k2  -  10k + 29)  + _ £ _ (k7  - 

256oc8(k  +  1)6 


30k6  +  241k5  -  629k4  -225k3  +  1836k2  +  2039k  +  623)  -  - 


512a9(K  +  l)9 


(3k8  -  91k7  +  635k6  -  946k5  -  1308k4  +  921k3  + 1  183k2  -  140k  -  257)  + 


(B.47) 


- £ - -(k9  -47k8  +  608k7  -2900k6 + 3328k5  +  10552k4  - 

1024a10(K  + 1)10 

10316k3  -34016k2  -  25109k  -  6101)  - H _ <3k10  -  144k9  + 

2048a11  (k  +  1)u 


1749k8  -  6484k7  +  744k6  +  20824k5  +  6536k4  -  18804k3  -  10539k2  + 


3584k  +  2531)  + _ £1 _ (k11  -  68k10  +  1305k9  -  9860k8  + 

4096a12(K  + 1)12 


27560k7  +  13328k6  -  149276k5  -73340k4  +  370775k3  +  569292k2  ♦ 


316995k +  63912), 
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in  -  i[_!+ _ l _ + _ t _ (k - 1) + _ t _ (k*-k-6)  + 

d  a  4a2(K  +  l)  8a3(K  +  l)2  16a4(K  +  l)3 


- "f. - (  k3  -  81c2  +  3k  +  4)  + _ ]_ - ( K*  +  6k3  -  14k2  +  54k  +  45)  + 

32a5(K  +  l)4  64a6(K  +  l)5 

t _ (k5  -  19k4  +  60k3  +  32k2  -  45k  -  29)  +  j. _ -(k6  - 

128a7(K  + 1)6  256a8(K  + 1)7 

15k6  +  8k4  +  271k3  -  256k2  -  844k  -  405)  + _ t _ -(k7  -  34k6 + 251k5  • 

512a9(K+  )fi 


309k4  -  803k3  +  22k2  +  615k  +  257)  + - £. _ -(ic8  -  28k7  +  124k6  + 

1024a10(K  + l)9 


(B.48) 


576k5  -  3008k4  -  1816k3  +9196K2  + 1  1860k  +  4023)  + _ tl - 

2048a11(K  +  l)10 


(k9  -  53k8 + 688k7  -  2540k6  -  836k5  ♦  9844k4  +  8596k3  -  4976k2  -  8193k  - 


2531)  - _ I _ (k10  -  45k9  +  430k8  +  150k7  -  12270k6  +  19534k5  + 

4096a12(K  +  l)u 


6507C  K4  -  44178k3  -  198071k2  -  161797k  -  42528)], 
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— 1  -  _ L  + _ 1 _ _ (k - 1)  + - — - (k2 -4k-  1) - 

d  k  +  1  a(K  +  l)2  2a2(K+l)3  4a3(K  +  l)4 


_ £ _ (k3-7k2  +  3k+3)  + _ £ _ (k4  -  12k3  +  20k2  +  24k  +  7)  - 

8a4(K  + 1)5  16a5(K  + 1)6 


£  -(k5  17k4  +  50k3  +  22k2  -  35k  -  21)  + _ £ _ -(k6  -  24k5  + 

32a6(K  +  l)7  64a7  (k  +  1)8 


1  19k4  -  36k3  -  301k2  -  244k  -  59)  - _ £ - -(k7  -  31k6  +  213k5  - 

128a8(K  +  l)9 


251k4  +  605k3  +  35k2  +  455k  + 183 )  + _ £ _ •(  k8  -  40k7  +  384k6  - 

256a9(K  +  l)10 


904k5  - 1  184k4  +  2320k3  +  4636k2  +  2768k  +  563)  - 1 - -(k®  - 

512a10(K  +  l)u 


49k8  +  596k7  -  2084k6  -  658k5  +  7458k4  +  6244k3  -  3796k2  -  5927k  - 1785) 


+ _ £1 _ (k10  -  60k9  +  935k8  -  4720k7  +  2840k6  +  23272k5  -  3  16k4 

1024au(K  +  l)12 


-  58720k3  -  70049k2  -  33180k  -  5795), 
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•  it — + _ l! _ - _ t — (k-d+ — "f. — 

d  k  +  1  4cx3(k  +  1)3  8a4(K  +  l)4  8a5(K  +  l)5 

_ ]_  (k3-6k2  +  k+4)  + _ I - (3k4  -32^  + 41k2  +  110k +  50)  - 

16a6(K  +  l)6  64a7(K  +  l)7 

^  -  -(3k5  -  41k4  +  99k3  +  99k2  -  86k  -  74)  + _ t _ (  k6  -  20k5 

128a8(K  +  l)8  64a9(K  +  l)9 

(B.50) 

+  86k4  +  30k3  -  306k2  -  364k  - 1 19)  - tl _ -(2k7  -  48k6  +  276k5  - 

256a10(K+l)10 

187k4  -  1006k3  -  170k2  +  760k  +  373)  + - tl - -(5k8  -  160k7  + 

1024a11(K  +  l)11 

1330k6  -  2340k5  -  6786k4  +  7168k3  +24634K2  +  19156k  +  4897)]. 

B.4  Leading  terms  of  the  integrands  in  Fredholm  kernels 

It  is  shown  in  Chapter  2  that  the  Fredholm  kernels  in  the  singular  integral 
equations  are  in  the  form  of  infinite  integrals  with  their  integrands  composed  of  two 
multiplicative  parts  (Eqns.’  (51)  through  (54)).  The  first  is  made  up  of  sums  of 
analytical  rational  functions  whose  numerators  and  denominators  are  determinants 
as  described  in  Chapter  2.  These  expressions  are  multiplied  by  either  a  sine  or  cosine 
function  that  makes  the  integrands  oscillating. 

In  Chapter  3,  we  choose  to  evaluate  each  infinite  integral  by  two  parts.  The 
first  part  is  an  integral  with  limits  of  integration  from  0  to  a  finite  value  "A"  (see 
Chapter  3),  and  the  integral  is  evaluated  numerically  by  straightforward  Gauss’ 
formula.  The  second  is  an  integral  with  limits  of  integration  from  "A"  to  infinity.  We 
have  chosen  to  evaluate  only  the  leading  terms  of  the  integrand  from  the  asymptotic 
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expansion,  together  with  the  oscillating  part,  in  closed  form.  The  remaining  part  is 
ignored,  provided  a  sufficiently  large " A "  is  chosen  so  that  contributions  to  the  integral 
from  terms  of  order  higher  than  what  could  be  extracted  asymptotically  were 
insignificant.  Convergence  of  the  evaluation  of  each  infinite  integral  is  achieved  by 
varying  the  value  "A"  and  by  noting  the  numerical  results  are  identical  to  a  certain 
digit. 

In  what  follows,  we  shall  give  the  leading  terms  of  asymptotic  expansion  of  the 
integrands  in  the  Fredholm  kernels  as  described  above.  They  are  derived  simply  from 
combining  terms  and  expressions  already  established  in  Eqns.  (B.47)  through  (B.50). 

Let  us  rewrite  Eqns.  (47)  through  (54)  as  follows, 


kn{x,t)  -  Jo-D1’1(a)sina(t-x)  da,  (B.51) 

kn(x ,t)  «  Z?j"2(ot)  cosa(£  -x)  da,  (B.52) 

k21(x,t)  *  Jq  D21(a)  cosa(t  -x)  da,  (B.53) 

k22(x,t)  *  J^D22(a)  sina(f  -x)  da,  (B.54) 


Dmn(  a) 

D;2(a) 


JSli[2a(hl 

2  D 


^31 

D 


K-l 


2a(x  + 1) 


-Mk+ix: 


'21 


D 


^41 

D 


K  + 1 


)], 


Jill  i  [2a(£li  +  EH)  -(k  +  wEn  -  ^li)], 

2  D  D  D  D 


(B.55) 


(B.56) 
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D2»  -  JS^.[2a(£ii-^!i)-(K-lK^l+£li)], 
2  D  D  D  D 


(B.57) 


DU  a)  -  -Jilii[2a(£“-f^-_L)-(K-lX^+£^-_L)]. 
2  D  D  2a  D  D  k  +  1 


(B.58) 


Let  du',  d12,  d2l'  and  d22  represent  the  leading  asymptotic  terms  of  Dn\  D12, 

d  d  d  d 

D21'  and  D22  respectively.  Substituting  — 11,  _1£,  —11  and  —11  fromEqns. 

d  d  d  d 

(B.47)  through  (B.50),  which  are  the  first  eleven  leading  asymptotic  terms  of  the 

respective  rational  expression  resulting  from  the  7  by  7  determinant  in  the  numerator 

and  the  8  by  8  determinant  in  the  denominator,  for  Pll,  Rll,  P~  and  — 22 

D  D  D  D 

in  Eqns.  (B.55)  through  (B.58),  we  obtain  a  set  of  truncated  polynomials  in  o’1.  Again, 

any  term  that  involves  _ li,  _ 11,  _ —  or  _ —  is  ignored  for  the  same 

D  D  D  D  * 


reason  given  earlier,  that  they  are  of  much  lower  order  terms  and  do  not  participate 
in  the  asymptotic  analysis  or  contribute  to  the  leading  terms.  The  expressions  for 
cTu(a),  cT i2(a),  cf2i(a)  and  cf22( a)  are  shown  in  the  following  few  pages. 
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d‘n{ a)  -  _L  + - £ - (k-1)(k-3)  + - £ - 

4a  8a2(K  +  l)2  16a3(K  +  l)3 


(k3  +  2k2  -  27k  - 12)  + _ £ _ (K  -  IXk3  -  8k2  +  17k  + 10)  + 

32a4(K  +  l)4 


_ £ _ (k5  -  5k4  -  60k3  +  160k2  +  235k  +  85)  + . £.  

64a5(K  +  l)6  128a6(K  + 1)6 


(k5  -  17k4  +  76k3  -  72k2  -  173k  -  71)  + . 


V7 


256a7(K  +  l)7 


-•(k7  -  16k3 


51k5  +  795k4  -  449k3  -  2926k2  -  2645k  -  741) 


512  a8(K  +  l)8 


(k  -  1)(k7  -30k6  +  241k5  -  629k4  -  225k3  ♦  1836k2  +  2039k  +  623)  + 


1024a9(K  +  l)9 


•(k9  -31k*  +«0k7  +  1860k6  -  7504k5  -8280k4  + 


24228k3  +  48824k2  +  3  1771k  +  7227)  + - £1 - -(k  - 1) 

2048a10(K  +  l)10 


(k9  -47k8  +  608k7  -2900k6  +  3328k5  +  10552k4  -  10316k3  - 


34016k2 -25109k -6101) 


Y11 


4096an(K  +  l)u 


.<k11-50k10  + 


445k8  +  2290k8  -31320k7  +  38896k®  +  195188k5  -  44732k4  - 


628313k3  -  775750k2  -394945k  -75502) 


n 


civ 


j.  i  o' 


(B.59) 
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d[2( a)  -  —  + - - - <k-1)  + - 1 - -(k2  -3k -8)  + 

4a  8a2(K  +  l)  16a3(K  +  l)2 


- £ - (k-1Kk2-5k-2)  + - £ - *(k4-10k3-2k2  + 

32a4(K  +  l)3  64a5(K  + 1)4 


82k  +  57)  + _ £ _ -(k4  -  14k3  +  26k2  +  38k  + 13)  + 

128a6(K  +  l)5 


_ £ _ -(k6  -  2  Ik6  +  66k4  +  253k3  -  558k2  - 1  164k  -  505)  + 

256a7(K  +  l)6 


_ i. _ -(k  -  1)(k®  -  27k6  +  148k4  -  45k3  -  452k2  -  404k  - 

512a8(K  +  l)7 


109)  + - £ - -(k8  -  36k7  +  276k6  +  48k6  -  3936k4  + 

1024a9  (k  +  1)8 


392k3  +  14556k2  +  15724k + 4975)  + - £1 - *(k  - 1) 

2048a10(K  +  l)9 


(K8  -  44k7  +  460k6  - 1  168k6  -  1648k4  +  3424k3  +  7316k2  +  4700k  + 


1039)  + _ £1 _ <k10  -  55k®  +  740k8  -  2190k7  -  10250k6  + 

4096au(K  +  l)10 


37786k®  +  64314k4 - 107782k3 -285651k2 -  209903k  -  52322) 


G?2"i(a)  -  -dx*2(a). 


(B.60) 


(B.61) 


175 


dl2(a.) 


-J—  + - — - (k  - 1)  + _ £ - *(k  -  4)  + 

4a  8a2(K+ 1)  16a3(K  +  l) 


Y4 

32a4(K  +  l)2 


(k-1)(k-6) 


_ £ _ -(k3  -  11k2  + 

64a5(K  + 1)3 


17k  +  33)  + _ t _ -(k  -  IXk3  -  15k2  +  33k  +  45)  + 

128a6(K  + 1)4 


- I _ -(k5  -  22k4  +  100k3  +  25k2  -  419k  -  305)  + 

256a7(K  + 1)5 


- 1 - *(k  -  IXk5  -  28k4  +  164k3  - 57k2  - 639k  - 

512a8(K  +  l)6 


405)  + - £ - -(k7  -  37k6  +  329k5  -601k4  + 

1024a9  (k  +  1)7 


1959k3  +  2831k2  +  6829k  +  3071)  + _ £1 _ -(k  - 1) 

2048a10(K  +  l)8 


<k7  +  45k6  -  489k5  -  1289k4  -  2  199k3  +  5279k2  +  9741k  + 


4023)  + - £1 - <k9  -  56k8  +  816k7  -  3646k6  - 

4096au(K  + 1)9 


1284k5  +  29710k4  +  7460k3  -  86570k2  -  100545k  -  32734) 


11  c' 
y>  c22 j 

j- 1  o' 


(B.62) 
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Appendix  C 

Formula  Derived  From  Sine  And  Cosine  Integrals 


/; 


sina(£-x)  ,  . _ ..  .  n 

-  da  «  sign  (t-x) 

o  a  2 


(C.7) 


Therefore 


si(A(t-x))  *  SitA(£-x))  -  signit  -x)—. 

2 


(C.8) 


Now  we  are  ready  to  derive  the  necessary  closed  form  expressions  for  integrals 
to  be  used  in  Chapter  3.  Integrating  Eqn.  (C.6)  by  parts  as 


f"  sina(*-x)  f"  — 

*'A  a  -'A  a2(t-x) 

and  after  simplifying,  we  obtain 


-  cosa(?-r)  ££a  m  _  cosa (t-x)  |“  _  cos A(t-x) 


a  (t-x)  U  A(t-x) 


(C.9) 


f"  cosaK~x)d a  .  .  (t-x)si(A(t-x)). 

a2  A 


(C.10) 


Again  integrating  Eqn.  (C.10)  by  parts  and  simplifying,  we  obtain 


J; 


sina(£  x\^a  m  L_iLcosA(£ -x)  +  sinA(f -x)  +  si(A(t -x)).  (C.ll) 

2A  2A2  2 


A  a3 


By  repeating  JA-.'  bove  process  severed  times,  we  can  deduce  the  following  identities, 


f“  C0Sa(t-x)rin  -  CQsA(t-x)T  ~X^'^2n  ~2^  + 

Ja  a2"  (2n-l)!A2n-2>*1 


(C.12) 


U  (2n  -  l)\A2n'2J  (2n-l)! 
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f-  da  -  coaAtt-x)E  -2j -2)1  + 

*/a  a2"'1  tr  (2n-2)!A2n-2^1 


(C.13) 


sinAtt  -x)£  1— *)2°  1)(2 n  2-/'_1)!  +  (-1)»^(< Z5^Hsi(A(t -x)). 
p  (2n  -2)!A2"'2-'  (2n-2)! 

where  7i  is  any  positive  integer.  Here  we  have  derived  expressions  for  infinite 
integrals  with  a  definite  lower  limit  of  integration  whose  integrand  involves  cosine 
functions  with  even  algebraic  decaying  power  functions  as  well  as  for  similar  integrals 
with  integrands  that  have  sine  functions  and  odd  algebraic  decaying  power  functions. 

Applying  the  above  procedure  to  the  cosine  integral  as  we  do  to  the  sine 
integral,  integrating  Eqn.  (C.5)  by  parts  and  simplifying,  we  obtain 


J;^ 


sina(f._x)^a  m  _Lsina«-x)  -  (t-x)Ci(A(t  -x)). 


(C.14) 


Once  again  repeating  the  integration  by  parts  several  times,  we  deduce  the  following 
identities, 


f“  -COSa{t^  da  -  cosA(f  -,)£  <-irHt-xro-'X2n-y-l)l  + 
Ja  a2"’1  P  (2n  -2)lA2n'2J 


sin ,  (-1  r<Lmia<AU-x». 


Pi  (2n  -2)\A2n'2J'1 


(2n  -2)! 


(C.15) 


-X) 


Ida  -  cosA(f-x)E  l:lrl(^^'1(2n-2^1)!  ♦ 
P  {2n  -1)LA2""2' 


sinAff-x)^^:1^:^^:^!  ♦  (-iriLfg^acAtf-x)). 
p  (2n  -  l)\A2n'2j'1  (2n  - 1)! 


(C.16) 
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again,  n  is  any  positive  integer.  Eqn.  (C.15)  and  (C.16)  show  expressions  for  infinite 
integrals  with  their  integrands  having  cosine  or  sine  functions  and  corresponding  odd 
and  even  algebraic  decaying  power  functions. 
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Appendix  D 

Approximation  Of  Step  Function  and  Logarithm  Function 


D.l  Gauss’  formula  and  its  approximation  of  a  step  function 

The  famous  Gauss’  formula  for  the  evaluation  of  integrals  of  bounded  functions 
is  as  follows  [23] 


£  m  dt  «  ^ 

i  «1 

where  t,  is  the  i-th  zero  of  the  Legendre  polynomial  Pn(t).  The  weight  function  w,(t) 
and  remainder  Rn  are  as  follows, 


w, 


(i  -th 


ip' 


Rn  -  _ 22^1("!)4 _ f2n)(Z,),  -1<$<1. 

(2n  +  l)[(2n)!]3 


(D.2) 


For  smooth  bounded  functions,  the  Gauss’  formula  provides  an  accurate  evaluation  of 
definite  integrals.  To  evaluate  certain  functions  that  have  discontinuities,  however, 
Gauss’  formula,  or  any  numerical  scheme  for  that  matter,  will  fail  if  proper  care  is  not 
taken  to  address  the  discontinuities.  In  the  case  of  a  step  function,  we  will  find  that 
straight  away  Gauss’  formula  fails  miserably.  The  reason  is  that  general  numerical 
schemes  which  use  polynomials,  and  in  the  case  of  Gauss’  formula  which  uses 
orthogonal  Legendre  polynomials,  are  simply  unable  to  properly  approximate  a  step 
function.  Compared  to  evaluating  integrals  with  a  step  function,  those  with  a 
logarithm  function  fare  much  better  when  using  the  Gauss’  formula.  But  the  error 
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is  still  significant.  We  shall  illustrate  this  by  using  the  Gauss’  formula  to  evaluate 
integrals  with  a  simple  step  function  and  functions  involving  step  function  and  simple 
power  functions  to  compare  them  with  the  values  computed  from  exact  analytical 
expression.  Integrals  involving  a  logarithm  function  and  a  power  function  will  be 
demonstrated  likewise. 

In  Chapter  3,  we  evaluated  infinite  integrals  with  a  definite  lower  limit  having 
integrands  with  power  functions  in  the  denominator  and  oscillatory  sine  or  cosine 
terms  of  the  form 

f-  da,  (-  cosa(t  -*>  da,  n  -  1  -  11.  (D.S) 

J  A  Q*  *)  A  GC*1 

From  the  closed  form  expressions  of  Eqn.  (D.3)  (see  Eqns.  (C.12)  ,(C.13),  (C.15)  and 
(C.16)  in  Appendix  C),  we  know  that  these  integrals  contain  terms  like  (t-x)  "  si(A(t-x)), 
and  (t-x)  "  Ci(A(t-x)).  From  Eqn.  (C.5)  in  Appendix  C,  we  see  that  Ci(A(t-x))  contains 
a  logarithm  term  like  log|£-x  |.  Furthermore,  from  Eqn.  (D.4),  we  know  that  si(A(t-x)) 
contains  a  step  function  sign(t-x).  This  implies  that  the  Fredholm  kernels  may  be 
expressed  as  follows 

kij(x,t)  «  jru-xr  [cij  n  sign(t  -x)  *  dtj  n  log|f-x|]  +  FtJ(tjc),  i,j  -  1,  2. 

n  »0 

where  ctJ  t  and  dLiJ  t  are  constants  and  Fi}(t,x)  are  smooth  bounded  functions.  The 
terms  involving  the  sign  function  in  the  summation  is  a  step  function  as  shown  in 
Figure  48.  The  n-1  term,  which  gives  (t-x)  sign(t-x),  is  continuous,  but  its  derivative 
is  not.  This  can  be  seen  from  the  presence  of  a  kink  in  Figure  49.  When  n  is  2  or 
larger,  each  term  becomes  increasingly  smoother,  even  though  there  is  a  step  function, 
and  it  becomes  easier  to  approximate  it  with  a  polynomial.  Figure  50  portrays  this 
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fact  by  the  plot  of  the  function  (t-xf  sign(t-x)  and  shows  that  even  with  the  presence 
of  a  step  function,  the  function  is  sufficiently  smooth  so  that  the  function  itself  and 
its  first  derivative  is  continuous,  only  its  second  derivative  is  discontinuous. 

Table  XXI  further  illustrates  the  fact  stated  above  by  the  tabulated  result  of 
numerical  integration  from  -1  to  1  of  functions  plotted  in  Figure  48  to  Figure  50  by 
three  different  schemes. 

(a) .  Use  a  20-point  Gauss’  formulas  over  the  limits  of  integration  from 

-1  to  1  disregarding  the  presence  of  step  function. 

(b) .  The  same  20-point  Gauss  rule,  but  splitting  it  into  two  integrals  at 
the  point  of  sign  change. 

(c) .  Exact  analytical  evaluation. 

We  make  several  observations  based  on  the  results  of  the  three  cases 
computed,  see  Table  XXI. 

(1) .  Completely  disregarding  the  presence  of  step  function  while 

integrating  numerically  could  make  the  integral  very  inaccurate, 
with  error  reaching  as  much  as  25%  in  the  cases  illustrated. 

(2) .  The  increasing  exponent  of  the  power  functions  serves  to  dampen 

the  effect  of  the  step  function  in  making  the  integrand  discontinuous. 
This  can  be  seen  in  the  increasing  digits  of  accuracy  for  the  integral 
as  the  power  goes  up. 

(3) .  Splitting  the  integral  in  two  when  there  is  a  step  function  in  the 

integrand  at  the  point  of  sign  change  is  perhaps  the  best  way  to 
evaluate  these  integrals.  This  can  be  seen  in  Table  XXI  that  the 
results  are  virtually  the  same  for  splitting  the  integral  and  the  exact 
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analytical  evaluation.  This  also  shows  the  superb  accuracy  Gauss’ 
formula  possesses  when  dealing  with  bounded  smooth  functions  in 
definite  integrals. 
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f(t)  =  sign(t-x). 
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(c).  Exact. 


D.2  Approximation  of  a  logarithm  function 

The  terms  in  Eqn.  (D.4)  that  contain  a  logarithm  function  do  not  behave  as 
badly  as  those  that  contain  a  step  function,  as  far  as  numerical  integration  is 
concerned.  As  a  matter  of  fact,  for  terms  like  (t-x)n  log|£-x|,  n  >  1,  the  function  is 
continuous  and  so  are  its  derivatives.  It  is  due  to  the  stronger  term  (t-xT  suppressing 
a  weak  singularity  at  (t-x)  =  0  of  a  logarithm  function.  This  can  be  seen  by  the 
behavior  of  the  functions  in  Figure  52  and  Figure  53.  Only  when  there  is  a  logarithm 
function  without  the  presence  of  power  function  do  we  see  a  singularity,  as  shown  in 
Figure  51. 

From  the  behavior  of  those  functions  that  contain  a  logarithm  term  as  shown 
in  Figure  51  through  Figure  53,  and  the  discussions  in  the  previous  section  regarding 
numerical  approximation  of  integrals,  we  could  very  well  anticipate  that  only  the 
logarithm  function  alone  would  need  particular  attention  in  numerical  computation. 
Those  integrals  that  have  terms  like  (t-x)n  log|*-x|,  n  >  1  pose  no  special  problem. 
Similar  to  what  has  been  done  in  the  previously  section,  a  table  is  compiled  to  list  the 
result  of  numerical  integration  from  -1  to  1  of  functions  plotted  in  Figure  51  to 
Figure  53,  again  using  three  different  schemes  for  comparison.  Indeed  we  find  such 
is  the  case  that  only  integration  of  the  log  term  itself  by  direct  application  of  Gauss’ 
formula  present  a  significant  error. 

Once  again,  we  make  the  following  observations  based  on  the  three  schemes 
used  for  these  functions  with  a  logarithm  term,  see  Table  XXII. 

(1).  The  inability,  or  difficulty  of  Gauss’  formula  to  appropriately  approximate 
functions  with  discontinuities  is  revealed  once  again  by  the  discrepancy  in 
the  result  of  direct  applying  the  formula  to  integrate  a  logarithm  function 
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and  the  exact  value.  Although  the  relative  error  in  evaluating  a  log 
function,  6%,  is  much  smaller  than  in  evaluating  a  step  function. 

(2) .  As  expected,  only  the  log  function  alone  causes  significant  error  in  the 
numerical  integration.  So  long  as  there  is  a  power  function  like  (f-x)  "  , 

n  >  1  to  suppress  the  logarithm  term,  Gauss’  formula  can  be  applied  directly 
without  significant  error. 

(3) .  Once  again,  splitting  the  integral  in  two  at  the  point  of  discontinuity  proves 
to  be  a  superior  way  of  handling  these  integrals. 

Unlike  the  case  of  a  step  function,  splitting  the  integral  when  integrating  a 
logarithm  function  from  -1  to  1  does  not  produce  as  good  an  accuracy.  This  is  because 
the  logarithm  function  is  unbounded  at  where  the  integral  is  being  separated. 
Artificially  separating  the  integral  at  point  of  singularity  forces  Gauss’  formula  to  take 
the  value  0  for  the  logarithm  instead  of  infinity.  But  the  error  introduced  is  still  small 
compared  to  if  Gauss’  formula  is  applied  over  the  whole  range  without  separation. 
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m 


Figure  51  log|f-x|,  x=0.3. 
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-0.4-1 


Table  XXII  Evaluation  of  /.j1  fU,x)  dt,  with  different  schemes  for  selected  x  values. 
f(t,x)  contains  a  logarithm  function. 


fU,x) 

Gauss 

Gauss/split  int. 

Exact 

(a) 

(b) 

(c) 

logtt-x) 

-1.807367 

-1.905593 

-1.908599 

(t-x)  log(t-x) 

-0.008847 

-0.009081 

-0.009083 

(t-x)2  log(t-x) 

-0.131007 

-0.130864 

-0.130864 

log(t-x) 

-1.825950 

-1.735370 

-1.738376 

x=0.5 

(t-x)  logft-x) 

-0.045311 

-0.042789 

-0.042792 

(t-x)2  log(t-x) 

0.038481 

0.038378 

0.038378 

log(t-x) 

-1.384621 

-1.456118 

-1.459124 

(t-x)  log(t-x) 

-1.121327 

-1.120933 

-1.120937 

(t-x)2  log(t-x) 

0.309208 

0.309268 

0.309268 

Note: 

(a) .  Gauss’  formula  with  20  Gauss  points. 

(b) .  Split  the  integral  in  two  at  the  point  (t-x)  =  0,  then  use 

Gauss’  formula  for  each  integral. 

(c) .  Exact. 
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Appendix  E 

Numerical  Evaluation  Of  Determinant 


E.l  Introduction 

In  Chapter  3,  we  solved  the  SIE  numerically  by  discretizing  it  at  the  collocation 
points.  At  each  collocation  point  we  evaluated  the  Fredholm  kernels  which  are  in  the 
form  of  infinite  integrals.  Those  integrals  are  each  separated  into  two  parts.  The  first 
part  is  obtained  by  using  closed  form  expressions.  The  second  part  is  a  definite 
integral  with  limits  of  integration  being  0  and  a  finite  number  "A",  which  we 
evaluated  by  using  Gauss’  formula. 

We  recall  in  Chapter  2  that  the  integrand  of  each  Fredholm  kernel  is  the 
algebraic  sum  of  four  terms  each  being  the  quotient  of  a  7  by  7  and  an  8  by  8 
determinant  (see  Eqns.  (51)  through  (54)).  From  the  discussions  in  Chapter  3  and 
Appendix  B,  we  realized  that  in  the  actual  numerical  computation  these  determinants 
are  to  be  evaluated  by  a  full  expansion.  Since  full  determinant  expansion  requires  a 
very  large  number  of  operations,  we  can  expect  this  is  where  the  major  computing 
effort  will  be  in  the  numerical  analysis. 

E.2  Floating  point  overflow 

In  order  to  learn  how  to  be  more  efficient  in  computing  the  8  by  8  and  7  by  7 
determinants,  let  us  write  again  the  8  by  8  determinant  shown  in  Eqn.  (B.3),  but  in 
a  slightly  different  format  as  shown  in  Eqn.  (E.l).  Following  discussions  in  Appendix 
B,  we  know  that  the  leading  terms  of  the  8  by  8  determinant  and  the  eight  7  by  7 
cofactors  appear  in  the  integrands  of  Fredholm  kernels  must  carry  a  term  like 
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e2a<fc‘ ***)  .  Depending  on  the  relative  value  of  a,  h1  and  h2,  this  e2a<h'  *  ^  term  can 
get  very  large  very  easily.  Floating  point  overflow  in  numerical  computations  is 
consequently  something  we  should  watch  out  for  in  this  determinant  expansion. 

Let’s  take  the  maximum  floating  point  number  prescribed  in  a  particular 
computer  to  be  R,  a  machine  constant.  This  number  is  hardware  dependent  as  well 
as  programming  language  dependent.  Let’s  further  restrict  the  programming 
language  used  to  be  FORTRAN,  which  is  the  programming  language  numerical 
procedures  in  this  study  is  written.  R  is  therefore  the  maximum  real  number,  single 
or  double  precision,  that  can  appear  anywhere  in  a  FORTRAN  program  before 
overflow  occurs.  In  the  8  by  8  determinant  that  we  are  dealing  with,  however,  all 
numerical  operations  are  in  complex  arithmetic.  The  maximum  floating  point  number 
that  can  appear  in  our  numerical  procedure  in  either  the  real  part  or  the  imaginary 

i 

part,  therefore,  becomes  R 7  .  The  numerical  computation  for  this  work  is  carried 
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out  on  a  Cyber  180  Model  850  where  the  R  value  allowed  in  FORTRAN  is  defined  to 
be  5.2xl01232  [32].  From  this  we  can  deduce  that  the  following  condition  must  exist 
or  a  floating  point  overflow  would  occur. 

e2 “*■**■)  £  Rt,  =>  a(h^h2)  Z  1  log R  -  709.61  (E-2) 

4 

Based  on  the  analysis  above,  we  now  know  that  in  evaluating  the  Fredholm 

kernels,  as  the  integration  variable  gets  larger,  there  is  increasing  danger  of  floating 

point  overflow.  One  way  to  avoid  the  floating  point  overflow  is  to  have  checkpoints 

in  the  numerical  procedure  so  that  whenever  a  particular  element  with  the 

exponential  term  gets  too  large,  the  entire  row  is  divided  by  a  constant  to  keep  it  in 

check.  This  division  on  rows  has  to  be  done  on  both  the  numerator  and  the 

denominator  so  that  the  quotient  is  unchanged.  But  as  this  division  is  carried  out, 

there  are  elements  on  the  same  rows  that  will  run  the  risk  of  floating  point  underflow. 

These  elements  are  those  that  carry  the  negative  exponential  as  shown  in  Eqn.  (E.l). 

To  remedy  the  situation  another  checkpoint  will  have  to  be  installed  to  set  those 

elements  under  risk  of  floating  point  underflow  to  zero. 

While  this  method  works,  there  is  another  way  of  computing  these 

determinants  more  efficiently.  We  shall  describe  it  as  follows.  The  value  —  log  R 

4 

in  Eqn.  (E.2)  serve  as  a  guide  below  which  full  expansion  should  be  used  to 
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numerically  compute  the  determinants.  When  a  (hx  +  h2)  >  —  log  R  ,  the  term 

4 

e2a (A,*v  ig  so  overwhelmingly  large  that,  as  a  result,  the  product  of  the  rest  of 

the  terms  of  the  fully  expanded  determinant  are  lost  in  the  trailing  digits  of  this  huge 

number  so  that  we  are  in  effect  computing  the  leading  terms  of  the  determinants.  If 

we  recognize  this  fact,  the  complicated  procedure  of  computing  all  the  terms  of  the 

determinant  shrunk  to  just  computing  the  leading  term  only,  that  is  the  term  with  the 

common  factor  e2a(h'  * Al)  .  When  this  happens,  we  don’t  compute  the  term 

e2a(h‘  *  h*1  in  the  numerator  nor  in  the  denominator  because  they  factor  out.  What 

follows  is  a  much  simplified  computation  of  the  determinants  exactly  as  described  in 

Appendix  B.  This  reflects  a  much  reduced  effort  in  our  numerical  computation. 

—  log  R  in  Eqn.  (E.2)  therefore  becomes  a  parameter  to  flag  between  full 
4 

expansion  and  leading  terms  expansion. 

To  avoid  any  possibility  of  floating  point  overflow  in  numerical  computation, 

we  purposely  keep  the  value  of  a-(ht  +  hj  smaller  than  —  log  R  by  some  margin. 

4 

Since  in  full  determinant  expansion,  we  have  eight  elements  to  contend  with  in  each 
product,  and  e2a(h'  *  h,)  constitutes  only  the  exponential  term  of  two  elements.  A 
limit  of  500.0  was  used  for  the  computation  in  this  work  instead  of  709.61  as  given  by 
Eqn.  (E.2). 

There  is  yet  another  interesting  point  one  may  make  as  a  result  of  this  floating 
point  overflow  discussion  made  above.  Take  the  integral  in  Fredholm  kernels  that  we 
have  chosen  to  ignore  in  Chapter  3,  namely 
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[Oi,(a)  -^(alJsinatt -x)  da, 
j~  [D2*2(a) -d  22(a)] sina(£  -x)  da, 

(E.3) 

[D1*2(a)-<f1*2(a)]cosa^-x)  da, 

[D2*i(oc)-^2*i(«)]cosa(^-x)  da. 

From  Chapter  2  we  know  that  D,  f's  are  the  algebraic  sum  of  Rl,  which  are 

quotients  of  full  expanded  determinants.  Likewise  d,  J  are  the  algebraic  sum  of 
d 

— ,  which  are  quotients  of  the  leading  terms  of  determinants.  From  the 
d 

discussions  stated  earlier,  we  conclude  that  if  the  integration  variable  a  gets  too  large, 
floating  point  overflow  can  occur.  Here  we  must  point  out  that  the  term  "large" 
depends  also  on  the  sum  (hj  +  hj.  For  large  (h}  +  hj,  a  can  be  quite  small  when 
floating  point  overflow  happens.  When  it  happens,  DtJ' s  are  evaluated  the  same  way 
as  dtJ  as  pointed  out  earlier,  that  is  by  leading  terms  expansion.  As  a  result  we  have 
no  way  to  reliably  evaluate  the  integrals  in  Eqn.  (E.3)  since  at  large  a,  we  can  not 
evaluate  DtJ'  exactly  due  to  limitation  of  the  machine  constant.  Putting  it  in  another 
way,  for  sufficiently  large  "A",  is  equal  to  d,'  as  far  as  numerical  computation  is 
concerned,  therefore  the  integral  is  zero.  This  proves  once  more  that  these  integrals 
can  be  ignored  as  we  have  done  in  Chapter  3. 

E.3  Case  of  thin  film  on  thick  substrate 

The  discussions  above  applies  to  all  cases  of  material  thicknesses  except  that 
which  we  are  going  to  address  in  this  section.  While  it  is  generally  true  that  when 
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the  order  of  magnitude  of  the  determinant  reaches  the  limit  of  machine  constant  (the 

i 

value  R  in  the  previous  section,  for  complex  arithmetic),  it  is  as  though  we  are 

computing  the  leading  terms  of  the  determinant  and  it  is  more  efficient  that  way.  But 

as  the  ratio  of  the  thickness  of  the  two  materials  —  gets  smaller,  it  reaches  a  point 

h  i 

where  the  numerical  scheme  of  computing  the  determinants  based  on  the  leading 
terms  expansion  no  longer  suffice.  Therefore  for  the  case  of  very  thin  film  on  thick 
substrate,  we  need  to  re-examine  the  numerical  procedure  because  the  usual  leading 
term  expansion  applicable  to  the  general  case  would  lead  to  difficulties  in  the 
convergence  of  solution  and  an  unsatisfactory  result.  We  shall  now  investigate  this 
situation  as  follows.  Consider  the  case  of  a  very  thin  nonhomogeneous  film  sis 
material  2 

h2  -  hhv  0E.4) 

Eqn.  (E.2)  becomes 


since  5  is  very  small.  We  now  consider  the  next  higher  order  term  in  the  determinant 
that  is  ignored  by  the  leading  term  expansion  but  a  term  that  becomes  increasingly 
important  as  —  becomes  smaller,  they  are  terms  with  a  factor  like  e 2ah'  (see  Eqn. 
(E.l)).  We  actually  do  not  know  under  what  circumstance  the  dropping  of  the  term  e2ah' 
will  affect  the  result.  Let  us  illustrate  by  assuming  conservatively  that  when  the 


second  highest  order  term  is  lsirger  than  1%  of  the  highest  order  term,  computing  the 
determinant  based  on  only  the  leading  terms  will  give  a  distorted  result.  In  other 
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words,  when 


e2aA>  >  1%  e2^'^, 


(E.6) 


the  solution  will  be  disturbed  if  terms  including  e2ah'  and  lower  are  dropped. 
Taking  logarithm  on  both  sides  of  the  inequality  and  simplifying,  we  obtain 


§  <  log  100 
2a/ix 


(E.7) 


As  mentioned  earlier,  as  e2aQ>'  * Al)  reaches  the  machine  constant,  we  resort  to 

leading  term  expansion  in  computing  the  determinants.  When  that  happens,  a  h1  is 

at  least  equal  to  JL  log  R  ,  see  Eqn.  (E.5).  Substituting  this  into  (E.7),  we  obtain 
4 

the  limiting  value  of  5  under  which,  based  on  the  assumption  we  have  made  in 

keeping  only  the  highest  order  term  when  the  second  highest  order  term  is  larger  than 

1%  of  the  highest  order  term,  using  the  leading  terms  expansion  to  compute  the 

determinants  is  no  longer  valid.  For  Cyber  850,  we  get  8  <  0.00324,  which  means  that 

at  approximately  —  >  300  ,  dropping  lower  order  terms  as  have  been  done  in 

k2 

leading  terms  expansion  will  affect  the  accuracy  of  our  solution. 

The  fix  for  these  borderline  cases  when  !ll  gets  large  is  to  simply  go  back 

h2 

to  the  less  efficient  full  determinant  expansion  and  employ  row- wise  division  to  avoid 

floating  point  overflow  as  well  as  setting  small  numbers  to  zero  to  avoid  floating  point 

underflow.  To  be  on  the  safe  side,  in  the  cases  that  we  have  computed  in  Chapter  4, 

h. 

two  worst  cases  where  —  equals  200  and  400  respectively  (hj  =  100  a,  h2  = 0.5  a, 

h2 

and  hj  =  100  a,  h2  = 0.25  a)  are  computed  using  full  determinant  expansion.  The  rest 


of  the  cases  were  computed  by  employing  the  leading  terms  expansion.  We  found  the 
results  to  be  quite  satisfactory  as  have  been  shown  in  Chapter  4. 


Appendix  F 

Behavior  of  Cauchy  Integral 


Let  there  be  a  function  of  position  <j y(t)  defined  on  an  arc  [a;,  a2].  ty(t)  is  said 
to  satisfy  the  Holder  condition  on  [av  a2],  if  for  any  two  points  tp  t2  on  the  arc 


£  A  |*2-*a  I”,  0<p^l, 


(F.l) 


where  A  is  a  positive  constant.  A  is  called  the  Holder  constant  and  ji  the  Holder 
index.  The  implication  that  the  function  $(t)  satisfies  the  Holder  condition  is  that  §(t) 
is  bounded  and  dty(t)/dt  is  integrable  within  the  closed  interval  [a  ,,  a2]. 

Consider  the  Cauchy  integral 

<b(z)  -  _L  f"  M.  dt,  (F.2) 

2x1  •'«>  t-z 

where  §(t)  has  at  most  an  integrable  singularity  and  may  be  expressed  as 

i?(t)  -  ML  +  -ML  +  <p- a),  o<i?e[Y1],  ^*3) 

(t -aj'  (t  -a2)r* 

<pM'(t),  k-1,2  satisfy  the  Holder  condition  on  the  closed  interval  of  [av  a2]  and  ^1’(a1) 
*  0,  b/faj)  *  0.  Substituting  Eqn.  (F.3)  into  Eqn.  (F.2),  <t >(z)  may  be  expressed  as 


<&(*)  -  J-  [« hiaj  f°’ - ± -  +  $l(a2)  f°’ 

2m  •*>.  (t-a^'it-z)  Ja' 


dt 


(A) 


(t-a2)yi(t-z) 

(I2) 


fa,  dt  +  pH  ^(O-Qafoa)  ^  +  j-a,  4>3«)  dt  ^ 

^  (t-djY'it-z)  h  ( t-a2)Ht-z )  J°> 

(i3)  a4)  a.) 


(F.4) 


Since,  satisfies  the  Holder  condition,  we  have 

|<t.;W-<D*(a)|  <  A  |f -a  I”,  A>0,  0<p<l.  (F.5) 

It  follows  that 

/,<A|_ L  f - * -  |<X| - 1 - 1,  (P.6) 

2*i  •**.  (t-aj'-'it-z)  (z-aj'-* 

where  A  and  K  are  both  positive  constants.  From  this  we  conclude  that  I3  is  of  higher 
order  than  ( z  -  a1)‘1'' .  By  a  similar  argument,  I4  is  of  higher  order  than  ( z  -  a2)~^  . 
§3( t)  also  satisfies  the  Holder  condition,  therefore  Is  has  at  most  a  log  singularity.  As 
for  Ij  and  4  we  shall  show  that  they  both  have  singularities  at  the  end  points  that 
is  stronger  than  any  of  I&  I4,  or  15.  Let  t0  be  a  point  on  (a,,  a2).  From  the  Plemelj 
formula  (see  [20])  we  obtain  the  boundary  value  of  7;  as  follows 


Ix(tor  -  Ii(toy  •  <1 >I(a1)(t0-al)T,I 

W  +  Ixu0y  -  f- - f lL - 

*i  {t-aj'it-tj 


(F.7) 


Define  the  function 
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on  any  definite  branch  that  is  single  valued  and  varies  continuously  on  the  interval 
(av  o2).  We  obtain  its  boundary  values  as 


n(t0r  -  <h(^--  (t0-air\ 

1  -e'2’“Tl 


Q(t°r  m  r 

1  -c'2**T' 


Qu0)*  - 


-Ti 


From  Eqn.  (F.7)  and  Eqn.  (F.9),  we  obtain 

-  tw-ac^o)]'  *  o. 

which  means  that  [I/z)  -  is  holomorphic  on  the  entire  plane. 
l,{z)  -  Q(z)  -  PJz), 


l,(z) 


4>i‘(ai) 

1  -e'2*n‘ 


{ z-a 1Yy'  +  P^z), 


<j>I(q1)e”Tl 

2 1  sin  7ryx 


(z  -di)"1' 


+  Pjiz). 


(F.9) 


(F.10) 

We  can  write 

(F.ll) 


where  P/z)  is  holomorphic.  To  obtain  the  boundary  value  of  lt(z)  on  [a,,  a2],  use  Eqn. 
(F.9)  and  the  Plemelj  formula 


w  -  -  W  ♦  Wl 


(F.12) 


_3_ ( t0-a j)^1  cotnyv 
2  i 


I2  can  be  obtained  in  a  similar  manner  in  terms  of  its  leading  term  and  a  holomorphic 
function.  By  Plemelj  formula  12  has  boundary  values  as  follows 


l2(t0y  -  /2(f0)"  -  Q>'2(a2)(t0-a2)r\ 

i2(t0r  +  i2(t0r  -  f* - - 

**  (t-a2)y'(t-t0) 


(F.13) 


Again,  define  a  function  'V(z)  that  is  single  valued  and  varies  continuously  on  any 
definite  branch 


*P(2)  - 


<j>2(a2) 


(1  -e2K,y,)(z  -a2)y' 


The  boundary  values  of  'V(z)  become 


vdor  -  .-fey  (t0-a2y\ 
l-e2*,y> 


1  _e2*‘YS 


T«0)*  -  ¥(t0r  -  ♦;(a2)(f0-atrk. 


(F.14) 


(F.15) 


From  Eqn.  (F.13)  and  (F.15),  we  conclude  that  I2(z)  -  'V(z)  is  holomorphic  and  write 
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I2(z)  -  'V(z)  -  P2(z), 


I2(z)  -  4>2-a.2).  (z  -a2r>  +  P2(z), 

l_e2-ra 


(F.16) 


<j>2(a2)e 


-**  Tj 


-te  -a,)'1'1  +  P,(2). 


2isin7T72 

Similar  to  Eqn.  (F.ll),  P/z)  is  a  holomorphic  function.  The  boundary  values  of  I/z) 
on  [a 2,  a2 ]  car  be  obtained  from  Eqn.  (F.15)  by  using  Plemelj  formula  again  as  follows 


/,««,)  -  4  tw +  w‘] 


(F.17) 


-L^jCaj)  (£0-a2HJ  cot7cy2. 
Zl 


We  have  now  shown  that  I2(z)  and  J/zj  have  higher  order  leading  terms  than  either 
I3(z)  or  I4(z).  As  for  Is(z)  which  has  a  log  singularity,  we  shall  assume  temporarily  that 
it  has  a  weaker  singularity  than  either  Ij(z)  or  I/z).  We  can  now  write  the  Cauchy 
integral  and  the  boundary  values  in  terms  of  its  leading  terms  as  follows 


<D(z) 


J_  f-  a, 

2ni  •*».  t-z 


213^717! 


( z  -a2V' 


2isinjr72 


(. z-a 2)~*  +  Piz). 


(F.18) 


-  _L<j>[(a1)cot7r71(^0 -aj)"’'1  -  _L<j>J(a2)cotjtY2(*0-a2)‘T’  +  P(t0).  (F.19) 

2  i  2i 

P(z)  is  bounded  everywhere  except  possibly  at  the  end  points  a,  b,  where  it  has 
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singularities  no  higher  than  those  of  the  leading  terms  in  Eqn.  (F.18). 

Multiplying  Eqn.  (F.19)  by  ( z  -  and  let  t„  — >  alt  we  obtain  the  following 

characteristic  equation 


^(ajlcotJCYj  _  n  (F.20) 

Ti  " 

Since  <^/(a,)  *  0,  we  have 


cotKYj  -  0,  Yi 


(F.21) 


where  Yj  represent  the  stress  singularity  at  end  a}.  We  obtain  Yi  "  —  as  the  only 

2 

acceptable  solution  since  0  <  Re  [y]  <  1,  At  the  end  a2  we  can  follow  a  similar 

procedure  giving  Yi  *  — ■  •  This  is  the  well-known  square  root  singularity  for  most 

2 

crack  problems.  It  also  justifies  the  assumption  we  made  earlier  that  I5  has  a  weaker 
singularity  than  that  of  It  and  I2,  since  log  singularity  is  weaker  than  square  root 
singularity.  From  the  discussions  we  had  so  far,  we  conclude  that  for  singular 
integral  equations  with  a  simple  Cauchy  kernel.  The  density  function  §(t)  may  be 


written  as 


4 >(*) 


r(t) 

(«-aj)1/2(a2-<)1/2 


<!>’(«)  -  ((ilttMaj,-*)1'2  +  i^Xf-aj)172  +  ^UXt -t)1*. 
where  ty'(t)  is  bounded  everywhere  on  [a,,  a2]. 


(F.22) 
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